优秀实验报告一.docx
- 文档编号:5418112
- 上传时间:2022-12-16
- 格式:DOCX
- 页数:35
- 大小:1.39MB
优秀实验报告一.docx
《优秀实验报告一.docx》由会员分享,可在线阅读,更多相关《优秀实验报告一.docx(35页珍藏版)》请在冰豆网上搜索。
优秀实验报告一
实验2:
π的计算
学院:
机械与动力工程学院姓名:
唐子彦学号:
515020910137
1、实验目的
通过求π的近似值,了解历史上计算π值的一些方法,包括刘徽割圆术、级数展开、数值积分和MonteCarlo法等.复习微积分中相关知识,比较它们的差异,了解计算方法对提高计算效率的意义.
2、实验内容
(1)利用鲁道夫割圆术计算π值
【问题】德国人鲁道夫一生计算圆周率,他同样是用圆的内接多边形逼近圆周.不过,他是从正方形开始成倍增加边数.试推导出他计算所采用的递推公式,然后求π的近似值到10位和20位有效数字.
【解】
图1.1
为通过“鲁道夫法”计算π的近似值,现编写M文件如下:
functionm2_8(m)
%该函数用于计算指定有效数字位数的π的近似值,m为有效数字位数
n=1;
PI=vpa(pi,100);
whileabs(PI-m2_7(n))>=10^(-m+1)
n=n+1;
end
disp(['迭代次数n=',num2str(n)])
vpa(m2_7(n),m)%仅显示π的近似值的指定位有效数字,下同
functionPI=m2_7(n)%该函数利用鲁道夫割圆术计算π的近似值,n为计算次数
a
(1)=sym(sqrt
(2));%为尽可能消除累计误差,此处采用符号运算
fori=1:
n-1
a(i+1)=sym(sqrt(2-sqrt(4-a(i)^2)));
end
S=sym(2^n*a(n));
PI=vpa(S,30);
为求出指定有效数字位数的π值,还需编写M文件如下:
运行结果如下:
图1.2
【答】由图1.2不难发现,“鲁道夫法”计算π值效率较高,n=16时即可得到π的10位有效数字,n=32时即可得到π的20位有效数字.然而,由于中途过程涉及开方等运算,因此计算较为复杂繁琐.
(2)利用幂级数展开式计算π值
【问题】简单公式
,Machin公式
,以及公式
.试验证上述三个公式(分别记为公式1、2、3),并利用反正切函数的幂级数展开式求π值,比较上述三公式的计算效率.此外,再找出一种利用幂级数展开式求π的方法并验证之.
【解】
为利用上述三个公式求出π值,现编写以下三个M文件:
functionPI=m2_10(n)
%该函数利用简单公式:
PI/4=arctan(1/2)+arctan(1/3)的幂级数展开式计算π
S=0;
fori=0:
n-1
S=S+(-1)^i/(2*i+1)*(1/(2^(2*i+1))+1/(3^(2*i+1)));
end
PI=vpa(4*S,50);
第一个M文件,对应于公式1:
functionPI=m2_11(n)
%该函数利用Machin公式:
PI/4=4*arctan(1/5)-arctan(1/239)的幂级数展开式计算π
S=0;
fori=0:
n-1
S=S+(-1)^i/(2*i+1)*(4*(1/5)^(2*i+1)-(1/239)^(2*i+1));
end
PI=vpa(4*S,50);
第二个M文件,对应于公式2:
functionPI=m2_12(n)
%该函数利用公式:
PI/4=arctan(1/2)+arctan(1/5)+arctan(1/8)的幂级数展开式计算π
S=0;
fori=0:
n-1
S=S+(-1)^i/(2*i+1)*((1/2)^(2*i+1)+(1/5)^(2*i+1)+(1/8)^(2*i+1));
end
PI=vpa(4*S,50);
第三个M文件,对应于公式3:
为比较它们计算指定有效数字位数的π值的效率,还需编写M文件如下(详见第五页):
functionm2_13(m)
%该函数用于计算对指定有效数字位数的π,利用不同幂级数展开式所需的项数
n1=1;n2=1;n3=1;
PI=vpa(pi,100);
whileabs(PI-m2_10(n1))>=10^-(m-1)
n1=n1+1;
end
whileabs(PI-m2_11(n2))>=10^-(m-1)
n2=n2+1;
end
whileabs(PI-m2_12(n3))>=10^-(m-1)
n3=n3+1;
end
disp(['有效数字位数m=',num2str(m)])
disp(['公式一所需的项数n1=',num2str(n1)])
disp(['公式二所需的项数n2=',num2str(n2)])
disp(['公式三所需的项数n3=',num2str(n3)])
disp('')
运行结果如下(每五个为一组):
图2.1图2.2
图2.3图2.4
从上述四幅图中可以看出,公式1、3的计算效率基本相同,而公式2的计算效率高于其他两个.
//利用arctan的Taylor级数展开式计算π的近似值
#include
#include
#definePI3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
constlongdoubleONE=1;
usingnamespacestd;
longdoublefun_1(intn);
longdoublefun_2(intn);
longdoublefun_3(intn);//下接第7页
然而,从有效数字位数m=15开始,所得结果中公式1、2、3所需项数不再增加,这可能是MATLAB本身的原因.个人猜测:
当MATLAB计算到一个与π极为相近的数时,可能将其自动补全为π,而没有继续计算.为试图解决该问题,可改用C++进行编程,所需的CPP文件如下:
intmain()//上接第6页
{
intn1=0,n2=0,n3=0;
intm;
cout<<"请输入有效数字位数m:
"< cin>>m; for(inti=1;i<=5;i++) { longdoubleepsilon=pow(ONE/10,m-1); while(abs(PI-fun_1(n1))>=epsilon) { n1++; } while(abs(PI-fun_2(n2))>=epsilon) { n2++; } while(abs(PI-fun_3(n3))>=epsilon) { n3++; } cout<<"为达到"< cout<<"π的近似值为: "< cout<<"为达到"< cout<<"π的近似值为: "< cout<<"为达到"< cout<<"π的近似值为: "< m++; cout< } return0; } longdoublefun_1(intn)//公式1利用π/4=arctan(1/2)+arctan(1/3) { longdoubleS=0; for(inti=0;i { longdoublea=pow(ONE/2,2*i+1); longdoubleb=pow(ONE/3,2*i+1); if(i%2==0) S=S+ONE/(2*i+1)*(a+b); else S=S-ONE/(2*i+1)*(a+b); }//下接第8页 return4*S;//上接第7页 } longdoublefun_2(intn)//公式2利用π/4=4arctan(1/5)-arctan(1/239) { longdoubleS=0; for(inti=0;i { longdoublea=pow(ONE/5,2*i+1); longdoubleb=pow(ONE/239,2*i+1); if(i%2==0) S=S+ONE/(2*i+1)*(4*a-b); else S=S-ONE/(2*i+1)*(4*a-b); } return4*S; } longdoublefun_3(intn)//公式3利用π/4=arctan(1/2)+arctan(1/5)+arctan(1/8) { longdoubleS=0; for(inti=0;i { longdoublea=pow(ONE/2,2*i+1); longdoubleb=pow(ONE/5,2*i+1); longdoublec=pow(ONE/8,2*i+1); if(i%2==0) S=S+ONE/(2*i+1)*(a+b+c); else S=S-ONE/(2*i+1)*(a+b+c); } return4*S; } 运行结果(每五个为一组)见第9页. 从运行结果来看,有效数字位数m=1~16均可得到正确的项数n1、n2、n3.然而当m=17时,或许是由于C++语言中longdouble类型的计算精度有限,无法进行高精度的浮点运算,导致循环条件恒为真,即程序进入死循环,无法得出正确结果.对于m>17的情形更是如此(参见图2.8). 图2.5图2.6 左图: 图2.7 下图: 图2.8 functionm2_28(m) %利用π/6=arcsin(1/2)及反正弦函数的Taylor级数计算π,m为有效数字位数 PI=vpa(pi,100); S=1/2;n=1;a=1;b=2; whileabs(PI/6-S)>=10^(-m+1) S=S+a/((2*n+1)*b)*(1/2)^(2*n+1); n=n+1; a=a*(2*n-1); b=b*(2*n); end disp(['为达到',num2str(m),'位有效数字,所需项数为n=',num2str(n)]) vpa(6*S,m) 为比较公式4与前三个公式的计算效率,现编写M文件如下: 输入图2.9所示命令行,运行结果如下: 图2.10 【答】对比以上四公式的计算结果不难发现,公式1、3计算效率大致相同且较低,公式2的计算效率最高,公式4的计算效率介于它们之间比公式1、3略高。 (3)数值积分计算π值 【问题】利用数值积分计算π,分别用“梯形法”和“Simpson法”精确到10位有效数字,再用“Simpson法”精确到15位有效数字. 【解】在本题中,“梯形法”和“Simpson法”的计算原理如下: 为利用数值积分计算π值,现编写如下两个M文件: functionPI=m2_17(n) %该函数利用“梯形法”进行数值计算,从而求出π的近似值,n为划分份数 x=0: 1/n: 1; y=1./(1+x.^2); S=2*sum(y)-1-0.5; PI=vpa(2*S/n,50); functionPI=m2_18(n) %该函数利用“Simpson法”进行数值计算,从而求出π的近似值,n为划分份数 x=0: 1/(2*n): 1; y=1./(1+x.^2); x1=1/n: 1/n: (n-1)/n; y1=1./(1+x1.^2); S1=sum(y1); S2=sum(y)-S1-1.5; S=1/(6*n)*(1.5+2*S1+4*S2); PI=vpa(4*S,50); 为计算指定有效数字位数的π值,还需编写以下两个M文件: (1)梯形法: functionm2_19(m) %该函数利用“梯形法”计算指定有效数字位数的π,以及最少的划分份数n n0=1; PI=vpa(pi,100); whileabs(PI-m2_17(n0))>=10^(-m+1) n0=2*n0; end a=n0/2;b=n0;n=(a+b)/2;%此处采用了“二分法”的思想,以便于快速找出n whileb-a>=10%本应为b-a>=1,此处是为了避免四舍五入产生的误差导致计算次数过多 ifabs(PI-m2_17(n))>=10^(-m+1) a=n; n=floor((a+b)/2); else b=n; n=floor((a+b)/2); end end whileabs(PI-m2_17(n))>=10^(-m+1) n=n+1; end disp(['划分份数至少为n=',num2str(n)]) vpa(m2_17(n),m) (2)Simpson法: functionm2_20(m) %该函数利用“Simpson法”计算指定有效数字位数的π,以及最少的划分份数n n=1; PI=vpa(pi,100); whileabs(PI-m2_18(n))>=10^(-m+1) n=n+1; end disp(['划分份数至少为n=',num2str(2*n)]) vpa(m2_18(n),m) 运行结果见第14页图3.1. 图3.1 【答】由图3.1知,利用“梯形法”和“Simpson法”计算10位有效数字的π值所需的最少划分份数分别为12910和20;利用“Simpson法”计算15位有效数字的π值所需的最少划分份数为92. 此外不难发现,当指定相同有效数字位数时,“Simpson法”所需的划分份数远少于“梯形法”.换言之,当划分份数一定时,“Simpson法”的计算精度远高于“梯形法”. 综上所述,“Simpson法”的计算效率远高于“梯形法”. (4)利用MonteCarlo法计算π值 【问题】 (1)试用MonteCarlo法计算π的5位有效数字; (2)设计方案试用计算机程序模拟Buffon实验. 【解】首先将MonteCarlo法和Buffon实验的原理进行一些简要说明. functionPI=m2_21(n)%MonteCarlo投点法 m=0; fork=1: n ifrand^2+rand^2<=1 m=m+1; end; end; PI=vpa(4*m/n,10); 由于取n=10000运行该程序1次和取n=1000运行该程序10次再求平均值的本质没有区别,故每次进行投点实验时,仅运行该程序一次.运行结果如下: 图4.1 从图4.1中可以看出,随着投点次数不断增加,所得的π值与真实值的误差总体上不断减小.本次实验中,当投点次数为10亿,可得π的5位有效数字. 然而,投点10万次,即可求出π的4位有效数字;而投点1000万次,仅可求出π的3位有效数字.这说明增加投点次数,不能说计算π的精度必定提高. 图4.2 图4.3 为利用计算机程序进行随机模拟,现编写M文件如下: functionm2_29(l,d,m)%该函数利用Buffon投针法计算π的近似值 %其中l为针长,d为平行线间距,m为试验次数 ifl/d<=0||l/d>1 disp('错误! 请重新输入l,d值。 ') else n=0; fori=1: m theta=pi*rand (1); h=-d/2+d*rand (1); ifabs(h)<=l/2*sin(theta) n=n+1; end end PI=vpa(2*l/d*m/n,10) end 运行结果如下: 图4.4 【答】尽管MonteCarlo法和Buffon实验的原理不复杂,操作也简单,但计算π效率很低,必须通过多次(往往达成百上千亿次)实验方可得到较精确的结果. (5)利用“Wallis公式”计算π值 【问题】试利用积分 推导公式: 并利用该公式计算π值,观察计算效率. 【解】 为利用Wallis公式计算π值,并观察其计算效率,现编写M文件如下: functionm2_23(m)%利用Wallis公式计算π的近似值,m为有效数字位数 S=1;n=0; PI=vpa(pi,100); whileabs(PI-2*S)>=10^(-m+1) n=n+1; ifmod(n,2)==1 a(n)=(n+1)/n; else a(n)=n/(n+1); end S=S*a(n); end disp(['为达到',num2str(m),'位有效数字,所需的项数n=',num2str(n)]) vpa(2*S,m) 运行结果如下: 图5.1 【答】图5.1表明,利用“Wallis公式”计算π值,每增加一位有效数字,所需的项数n就要提高一个数量级.因此,“Wallis公式”计算π的效率较低,大致和反正切函数公式 的幂级数展开式的计算效率处于同一数量级. (6)e的计算 【问题】试利用多种方法计算自然对数的底数e的近似值. 【解】以下四种方法均可用于计算e值. (1)利用公式1: . functionm2_24(m)%利用(1+1/n)^n计算e的近似值E,m为有效数字位数 E=exp (1); n=1; whileabs(E-(1+1/n)^n)>=10^(-m+1) n=n+1; end disp(['为达到',num2str(m),'位有效数字,所需的n值为',num2str(n)]) e=vpa((1+1/n)^n,m) 证: 这是e的定义式,故证明从略.现编写M文件如下: 输入以下命令: 图6.1.1 运行结果如下: 图6.1.2 (2)利用公式2: . 现编写M文件如下: functionm2_27(m)%利用(1+1/n)^(n+1)计算e的近似值E,m为有效数字位数 E=exp (1); n=1; whileabs(E-(1+1/n)^(n+1))>=10^(-m+1) n=n+1; end disp(['为达到',num2str(m),'位有效数字,所需的n值为',num2str(n)]) e=vpa((1+1/n)^(n+1),m) 输入以下命令: 上图: 图6.2.1 运行结果如下: 下图: 图6.2.2 (3)利用公式3: . functionm2_25(m)%利用Taylor级数计算e的近似值,m为有效数字位数 E=vpa(exp (1),100); e=0;n=0; whileabs(E-e)>=10^(-m+1) e=e+1/factorial(n); n=n+1; end disp(['为达到',num2str(m),'位有效数字,所需Taylor级数项数n=',num2str(n)]) e=vpa(e,m) 现编写M文件如下: 输入以下命令: 上图: 图6.3.1 运行结果如下: 下图: 图6.3.2 (4)利用公式4: 现编写M文件如下: functionm2_26(m,k) %利用(1+1/n+k/n^2)^n计算e的近似值E,m为有效数字位数,k为给定的正常数 E=exp (1); n=1; whileabs(E-(1+1/n+k/n^2)^n)>=10^(-m+1) n=n+1; end disp(['为达到',num2str(m),'位有效数字,所需的n值为',num2str(n)]) e=vpa((1+1/n+k/n^2)^n,m) 输入以下命令: 上图: 图6.4.1 所得结果如下: 下图: 图6.4.2 另外,为探究参数k对于计算效率的影响,进行了下述实验,运行结果如下: 图6.4.3 【答】对比图6.1.2、图6.2.2、图6.3.2、图6.4.2,不难发现公式3的计算效率最高,每提高一位精度只需增加一到两项,且随着n! 的爆炸增长,该公式计算e的效率会越来越高;公式1、2、4的效率彼此接近且较低,若要求的有效数字位数较多,则公式4的计算效率更高. 观察图6.4.3可以发现: 对于公式4,当m较小时,参数k越大,计算效率越低(对于较大的m值,由于计算机性能等原因,难以进行实验). 3、实验心得 通过本次实验,我对于π的几类常见计算方法及其计算效率有了初步的了解,对于MATLAB等编程软件的使用也不再陌生,同时也复习回顾了所学的微积分及概率等知识,感觉收获颇丰. 然而,在本次实验中,我也遭遇了一个问题,即所谓的“MATLAB自动补全功能”(参见第6页中部).为此,我做了两个小实验进行检验. 小实验1: 利用 的级数展开计算π,取n=20.编写M文件如下: functionS=m2_14(n) S=0; fori=1: n ifmod(i,2)==0 S=S-1/(2*i-1)*(1/(2^(2*i-1))+1/(3^(2*i-1))); else S=S+1/(2*i-1)*(1/(2^(2*i-1))+1/(3^(2*i-1))); end end S=vpa(4*S,30);%观察30位有效数 注: 上述M文件与《数学实验》(第二版)第22页下部的程序完全相同. 运行结果如下: 奇怪的是,相同的程序却得出了不同的结果.更令人意外的是,ans与π的前30位竟完全吻合,我认为这可能是MATLAB升级后新增了“自动补全功能”. 小实验2: functionPI=m2_15(n)%利用Ramanujan公式计算π的近似值 S=0; fori=0: n S=S+factorial(4*i)*(1103+26390*i)/((factorial(i))^4*396^(4*i)); end S=2*sqrt (2)/9801*S; PI=vpa(1/S,100); 利用拉马努金公式 计算π值,编写M文件: 运行结果如下: 令人震惊的是,借助拉马努金公式,仅用n=0和n=1两项即得到了π的100位有效数字,这绝不可能.因此我猜测,这应该是MATLAB的“自动补全功能”. 不论所谓的“自动补全功能”是否存在,这一现象直接导致在第二部分实验中,当要求的有效数字位数m较大时,利用MATLAB无法得出正确的所需项数n.
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 优秀 实验 报告