- 1、本文档共22页,可阅读全部内容。
- 2、原创力文档(book118)网站文档一经付费(服务费),不意味着购买了该文档的版权,仅供个人/单位学习、研究之用,不得用于商业用途,未经授权,严禁复制、发行、汇编、翻译或者网络传播等,侵权必究。
- 3、本站所有内容均由合作方或网友上传,本站不对文档的完整性、权威性及其观点立场正确性做任何保证或承诺!文档内容仅供研究参考,付费前请自行鉴别。如您付费,意味着您自己接受本站规则且自行承担风险,本站不退款、不进行额外附加服务;查看《如何避免下载的几个坑》。如果您已付费下载过本站文档,您可以点击 这里二次下载。
- 4、如文档侵犯商业秘密、侵犯著作权、侵犯人身权等,请点击“版权申诉”(推荐),也可以打举报电话:400-050-0827(电话支持时间:9:00-18:30)。
- 5、该文档为VIP文档,如果想要下载,成为VIP会员后,下载免费。
- 6、成为VIP后,下载本文档将扣除1次下载权益。下载后,不支持退款、换文档。如有疑问请联系我们。
- 7、成为VIP后,您将拥有八大权益,权益包括:VIP文档下载权益、阅读免打扰、文档格式转换、高级专利检索、专属身份标志、高级客服、多端互通、版权登记。
- 8、VIP文档为合作方或网友上传,每下载1次, 网站将根据用户上传文档的质量评分、类型等,对文档贡献者给予高额补贴、流量扶持。如果你也想贡献VIP文档。上传文档
查看更多
统计计算统计模拟π的估计.
统计模拟作业
班级 学号 姓名
作业一:π的估计
1. 分析方法
从十七世纪中叶起,人们开始用更先进的分析方法来求∏的近似值,其中应用的主要工具是收敛的无穷乘积和无穷级数,下面列举一些用此类方法求∏近似值的实例。
(1)由公式
推出=4
编写程序:
syms k
x=symsum((-1)^k/(2*k+1),k,0,10)
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.146567747182986
…
K π 10 3.232315809405593 50 3.161198612987050 100 3.151493401070990 200 3.146567747182986
(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.067703806643498
增加k的值
K=20,π≈3.103516961539230
K=50,π≈3.126078900215409
K=100,π≈3.133787490628159
K=10000,π≈3.141514118681864
K=1000000,π≈3.141591868191880
…
k π 10 3.067703806643498 20 3.103516961539230 50 3.126078900215409 100 3.133787490628159 10000 3.141514118681864 1000000 3.141591868191880
(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.14000000000000000000000000000
泰勒级数法
函数的泰勒级数展开式为:
将代入上式有
.
利用这个式子就可以求出的值了。
编写程序:
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.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);
end
文档评论(0)