统计计算统计模拟π的估计

1、统计模拟作业班级学号姓名作业一：的估计1. 分析方法 从十七世纪中叶起，人们开始用更先进的分析方法来求的近似值，其中应用的主要工具是收敛的无穷乘积和无穷级数，下面列举一些用此类方法求近似值的实例。(1)由公式 推出=4编写程序：syms kx=symsum(-1)k/(2*k+1),k,0,10)y=4*x得出当k=10时，3.232315809405593编写程序syms kx=symsum(-1)k/(2*k+1),k,0,20)y=4*x得出当k=20时，3.189184782277595依次，加大k的值K=50，3.161198612987050K=100，3.151493401070990K=200，3.146567747182986K103.232315809405593503.1611986129870501003.1514934010709902003.146567747182986（2）沃里斯(Wallis)方法在积分学中我们经常会遇到如下的沃利斯(Wallis)公式。沃利斯(Wallis)公式揭示了与整数之间的一种很不寻常的关系。但在实际学习中很少注意到沃利斯(Wa

2、llis)公式,更不会关注它的应用。实际上,沃利斯(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);endy=2*x得k=10时，3.067703806643498增加k的值K=20，3.103516961539230K=50，3.126078900215409K=100，3.133787490628159K=10000，3.141514118681864K=1000000，3.141591868191880k103.067703806643498203.103516961539230503.1260789002154091003.133787490628159100003.14151411868186410000003.14

3、1591868191880（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; endendp=vpa(4*m/n,30)程序运行结果：p =3.14000000000000000000000000000（4） 泰勒级数法函数的泰勒级数展开式为：将代入上式有.利用这个式子就可以求出的值了。编写程序：i=1;n=1000;s=0;for i=1:n s=s+(-1)(i-1)/(2*i-1);endp=vpa(4*s,3

4、0)程序运行结果：p =3.14059265383979413499559996126当取的值为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);ends=s1+s2;p=vpa(4*s,30)程序运行结果：p =3.14159265358979323846264338328 显然，级数收敛越快，取同样的值可以得到更高的精度。以同样的方法，能得出，程序和上面的一样。这样的近似值可以精确到几百位。（6） 麦琴公式由推出=4()，这个公式由英国天文学教授约翰麦琴于1706年发现。他利用这个公式计算到了100位的圆周率。麦琴公式每计算一项可以得到1.4位的十进制精度。因为它的计算过程中被乘数

5、和被除数都不大于长整数，所以可以很容易地在计算机上编程实现。编写程序：syms n;f1=(-1)(n-1)*(1/5)(2*n-1)/(2*n-1);f2=(-1)(n-1)*(1/239)(2*n-1)/(2*n-1);ans1=symsum(f1,n,1,28);ans2=symsum(f2,n,1,28);ans=vpa(4*(4*ans1-ans2),100)得3.141592653589793238462643383279502884197130451046268578972203255663716036677133432949011735665451127（7）用下列公式 求值，并了解的其他公式。编写程序：digits(10000);N=100000000;p=1;for m=1:N;m1=m*m;m2=2*m;n=(4*m1)/(m2-1)*(m2+1);p=p*n;endp=vpa(p*2);改变N的值，分析得到的不同值的大小与N的关系。PI=3. 141592653589793238462643383279502884197169399375105820974944

6、592307816406286208998628034825 当N=10*104时，=3.14151411868195662435709891724400222301483154296875当N=10*108时，=3.141592643066262180440162410377524793148040771484375 在使用MATLAB的过程中发现，当N值越大的时候，PI的值精确度越高，然而N大到一定范围后得出结果的过程也就变得非常慢，并且在10*108后每增加10倍后PI没有很大的变化，并且运算速度也变得非常慢了。2.概率方法 取一个二维数组 (x,y) ,取一个充分大的正整数n,重复n次，每次独立的从（0,1）中随机地取一对数x和y，分别检验x的平方和y的平方之和小于等于1是否成立。设n次试验中等式成立的共有m次，令4mn。但这种方法很难得到的较好的近似值。编写程序：m=0;for i=1:100000x=rand;y=rand;if x2+y2=1;m=m+1;else end end4*m/100000得3.136000000000000n=100000时，3.1399200

7、000000003. 数值积分方法 半径为1的圆的面积是。以圆心为原点建立直角坐标系，则圆在第一象限的扇形是由与轴，轴所围成的图形，扇形的面积。只要求出扇形的面积，就可得出的值。而扇形面积可近似等于定积分的值。 对于定积分的值，可以看做成曲线与轴，所围的曲边梯形的面积。把分成等分，既得个点，组成个小区间，每一个小区间与轴，所围成的图形是一个小曲边梯形。而梯形的面积计算公式是，对于第个小曲边梯形有上底为，下底为。所有小梯形的高都为。所以第个小曲边梯形的面积为。曲边梯形的总面积即定积分的值就是所有小梯形的面积总和。 为了避免根号，我们也可以利用积分得出的值。我们可以利用对求曲边梯形的面积来得出定积分的值，从而得出的值。设分点x1,x2,xn-1将积分区间0,1分成n等分。所有的曲边梯形的宽度都是h=1/n。记yi=f(xi).则第i个曲边梯形的面积A近似地等于梯形面积，即：A=(y(i-1)+yi)h/2。将所有这些梯形面积加起来就得到：A2/n2(y1+y2+yn-1)+y0+yn编写程序：n=10000;x=linspace(0,1,n+1);y=1./(1+x.2);h=1/n;an

8、s=vpa(4*h*trapz(y),11)得3.1415926519n=100000时，编写程序：n=100000;x=linspace(0,1,n+1);y=1./(1+x.2);h=1/n;ans=vpa(4*h*trapz(y),11)得3.1415926536也可利用积分公式编写程序：n=10000;x=linspace(0,1,n+1);y=sqrt(1-x.2);h=1/n;ans=vpa(4*h*trapz(y),11)得3.1415914776n=100000时，得3.1415926164结果分析对于分析方法来说，一般而言逼近速度太慢，运算庞大，对速度造成了很大影响。相对来说麦琴和泰勒级数法逼近的速度大大增加，得到的近似值精度较高。对于概率方法来说，随机性很大，同一个实验次数，得出的并不相同，有时差别还会很大，所以这种方法很难得到较好的近似值。 对于数值积分法来说，计算量较大，运算慢，不如分析方法中的麦琴和泰勒级数法好，泰勒级数法的精确度较高。作业二：正态分布随机数的产生方法一：平均值为a，标准偏差为b，长度为N的正态分布MATLAB程序如下：m=input(请输入平均值：);n=input(请输入标准差：);t=input(请输入数据长度：);%产生正态分布的随机数for i=1:t a=rand; b=rand; X(i)=sqrt(-2)*log(a)*cos(2*pi*b); Y=X*n+m;enddisp(Y);%求平均值和标准差M=mean(Y);N=std(Y);disp(M);disp(N);%将数据写入文本文件fid=fopen(noise.dat,w);Z=Y;fprintf(fid,%ft,Z);fclose(fid);%绘图histfit(Y);xlabel(随

《统计计算统计模拟π的估计》由会员枫**分享，可在线阅读，更多相关《统计计算统计模拟π的估计》请在金锄头文库上搜索。