计算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和阶数nr=taylor(f,n,v,a),指定自变量v、阶数n，计算f在a的级数(3)级数法由于利用arctanx的幂级数展开法的收敛较慢，可采用公式的计算来求pi值。(4)特殊公式（BBP）四、设计算法、编程求解4.1数值积分法梯形公式Matlab代码：formatlongx=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=103.104518326248318n=1003.140417031779045n=5003.141487477002141n=10003.141555466911028n=100003.1415914776113244.2幂级数法Matlab代码：(1)formatlongsymsxf=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=103.339682539682540n=1003.121594652591011n=5003.137592669589472n=10003.139592655589785n=100003.141392653591792(2)formatlongsymsxf=atan(x);t=taylor(f,10,x,0);x=1/5;s1=eval(t);x=1/239;s2=eval(t);pi=16*s1-4*s2当n=10时，pi=3.141592682404399级数法N=10的pi值3.3396825396825403.1415926824043994.3迭代法Matlab代码：formatlonga=1;b=1/sqrt(2);s=1/2;forn=1:1:10n,y=a;a=(a+b)/2;b=sqrt(b*y);c=a^2-b^2;s=s-2^n*c;pi=2*a^2/send结果：迭代法Pi值n=13.187672642712109n=23.141680293297657n=33.141592653895460n=43.141592653589831n=53.141592653589880n=63.141592653589978n=73.141592653590173n=83.141592653590564n=93.141592653591346n=103.1415926535929094.4蒙特卡罗方法Matlab代码：formatlongs=0;n=10;%n=100;n=1000;n=10000;n=100000;n=1000000fori=1:na=rand(1,2);ifa(1)^2+a(2)^2=1s=s+1;endendpi=4*s/n结果：随机模拟法Pi值n=102n=1003.240000000000000n=10003.132000000000000n=100003.135600000000000n=1000003.138480000000000n=10000003.1400240000000004.5BBP公式formatlongsymsxy=1/16^x*(4/(8*x+1)-2/(8*x+4)-1/(8*x+5)-1/(8*x+6));s=0;forx=0:1:10s1=eval(y);x,s=s+s1end结果：BBP公式Pi值n=03.133333333333333n=13.141422466422466n=23.141587390346582n=33.141592457567436n=43.141592645460337n=53.141592653228088n=63.141592653572881n=73.141592653588973n=83.141592653589752n=93.141592653589791n=103.141592653589793五、分析求解结果由上表可知，蒙特卡罗方法计算出的pi值与真实值的误差相差较大并且收敛速度很慢；对于级数法，但由于所选择的的级数方法、公式不同，得到的结果也就不同，收敛速度较慢，而的收敛速度就较快；数值积分法和迭代法准确度较高，但数值积分法的收敛速度没有迭代法快、精度高，所以一般情况下采用迭代法求近似值较准确。对于特殊BBP公式，其打破了传统的计算方法，可以直接计算pi的任何第n位数，而不是先计算前面的n-1位数。随着学习研究，利用特殊的公式计算明显地提高了pi的计算值的精确度。六、参考文献[1]刘慧颖.MTALABR2007基础教程[M].北京：清华大学出版社，2008.7,200[2]李庆扬，王能超，易大义.数值分析[M].北京：清华大学出版社，2008.7，187-190[3]数学实验课件