优秀实验报告一.docx

优秀实验报告一.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.