优秀实验报告一.docx

1、优秀实验报告一实验2：的计算学院：机械与动力工程学院 姓名：唐子彦 学号：5150209101371、实验目的通过求的近似值，了解历史上计算值的一些方法，包括刘徽割圆术、级数展开、数值积分和Monte Carlo法等. 复习微积分中相关知识，比较它们的差异，了解计算方法对提高计算效率的意义.2、实验内容（1）利用鲁道夫割圆术计算值【问题】德国人鲁道夫一生计算圆周率，他同样是用圆的内接多边形逼近圆周. 不过，他是从正方形开始成倍增加边数. 试推导出他计算所采用的递推公式，然后求的近似值到10位和20位有效数字.【解】图1.1为通过“鲁道夫法”计算的近似值，现编写M文件如下：function m2

2、_8(m) %该函数用于计算指定有效数字位数的的近似值，m为有效数字位数n=1;PI=vpa(pi,100);while abs(PI-m2_7(n)=10(-m+1) n=n+1;enddisp(迭代次数n=,num2str(n)vpa(m2_7(n),m) %仅显示的近似值的指定位有效数字，下同function PI=m2_7(n) %该函数利用鲁道夫割圆术计算的近似值，n为计算次数a(1)=sym(sqrt(2); %为尽可能消除累计误差，此处采用符号运算for i=1:n-1a(i+1)=sym(sqrt(2-sqrt(4-a(i)2);endS=sym(2n*a(n);PI=vpa(

3、S,30);为求出指定有效数字位数的值，还需编写M文件如下：运行结果如下： 图1.2【答】由图1.2不难发现，“鲁道夫法”计算值效率较高，n=16时即可得到的10位有效数字，n=32时即可得到的20位有效数字. 然而，由于中途过程涉及开方等运算，因此计算较为复杂繁琐.（2）利用幂级数展开式计算值【问题】简单公式，Machin公式，以及公式. 试验证上述三个公式（分别记为公式1、2、3），并利用反正切函数的幂级数展开式求值，比较上述三公式的计算效率. 此外，再找出一种利用幂级数展开式求的方法并验证之.【解】 为利用上述三个公式求出值，现编写以下三个M文件：function PI=m2_10(n)

4、%该函数利用简单公式：PI/4=arctan(1/2)+arctan(1/3)的幂级数展开式计算S=0;for i=0:n-1 S=S+(-1)i/(2*i+1)*(1/(2(2*i+1)+1/(3(2*i+1);endPI=vpa(4*S,50);第一个M文件，对应于公式1：function PI=m2_11(n) %该函数利用Machin公式：PI/4=4*arctan(1/5)-arctan(1/239)的幂级数展开式计算S=0;for i=0:n-1 S=S+(-1)i/(2*i+1)*(4*(1/5)(2*i+1)-(1/239)(2*i+1);endPI=vpa(4*S,50);第

5、二个M文件，对应于公式2：function PI=m2_12(n) %该函数利用公式：PI/4=arctan(1/2)+arctan(1/5)+arctan(1/8)的幂级数展开式计算S=0;for i=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);endPI=vpa(4*S,50);第三个M文件，对应于公式3：为比较它们计算指定有效数字位数的值的效率，还需编写M文件如下（详见第五页）： function m2_13(m)%该函数用于计算对指定有效数字位数的，利用不同幂级数展开式所需的项数n1=1;n2=1;n3=

6、1;PI=vpa(pi,100);while abs(PI-m2_10(n1)=10-(m-1)n1=n1+1;endwhile abs(PI-m2_11(n2)=10-(m-1)n2=n2+1;endwhile abs(PI-m2_12(n3)=10-(m-1)n3=n3+1;enddisp(有效数字位数m=,num2str(m)disp(公式一所需的项数n1=,num2str(n1)disp(公式二所需的项数n2=,num2str(n2)disp(公式三所需的项数n3=,num2str(n3)disp( )运行结果如下（每五个为一组）： 图2.1 图2.2 图2.3 图2.4 从上述四幅图

7、中可以看出，公式1、3的计算效率基本相同，而公式2的计算效率高于其他两个. /利用arctan的Taylor级数展开式计算的近似值#include #include #define PI 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068const long double ONE=1;using namespace std;long double fun_1(int n);long double fun_2(int n);long double

8、 fun_3(int n); /下接第7页然而，从有效数字位数m=15开始，所得结果中公式1、2、3所需项数不再增加，这可能是MATLAB本身的原因. 个人猜测：当MATLAB计算到一个与极为相近的数时，可能将其自动补全为，而没有继续计算. 为试图解决该问题，可改用C+进行编程，所需的CPP文件如下：int main() /上接第6页 int n1=0,n2=0,n3=0; int m; cout请输入有效数字位数m：m; for(int i=1;i=epsilon) n1+; while(abs(PI-fun_2(n2)=epsilon) n2+; while(abs(PI-fun_3(n3

9、)=epsilon) n3+; cout为达到m位有效数字，公式1所需的项数n1=n1endl; cout的近似值为：setprecision(m)fun_1(n1)endl; cout为达到m位有效数字，公式2所需的项数n2=n2endl; cout的近似值为：setprecision(m)fun_2(n2)endl; cout为达到m位有效数字，公式3所需的项数n3=n3endl; cout的近似值为：setprecision(m)fun_3(n3)endl; m+; coutendl; return 0;long double fun_1(int n) /公式1利用/4=arctan(1

10、/2)+arctan(1/3) long double S=0; for(int i=0;in;i+) long double a=pow(ONE/2,2*i+1); long double b=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页 return 4*S; /上接第7页long double fun_2(int n) /公式2利用/4=4arctan(1/5)-arctan(1/239) long double S=0; for(int i=0;in;i+) l

11、ong double a=pow(ONE/5,2*i+1); long double b=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); return 4*S;long double fun_3(int n) /公式3利用/4=arctan(1/2)+arctan(1/5)+arctan(1/8) long double S=0; for(int i=0;i17的情形更是如此（参见图2.8）. 图2.5 图2.6 左图：图2.7 下图：图2.8function m2_28(

12、m)%利用/6=arcsin(1/2)及反正弦函数的Taylor级数计算，m为有效数字位数PI=vpa(pi,100);S=1/2;n=1;a=1;b=2;while abs(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);enddisp(为达到,num2str(m),位有效数字，所需项数为n=,num2str(n)vpa(6*S,m) 为比较公式4与前三个公式的计算效率，现编写M文件如下： 输入图2.9所示命令行，运行结果如下： 图2.10【答】对比以上四公式的计算结果不难发现，公式1

13、、3计算效率大致相同且较低，公式2的计算效率最高，公式4的计算效率介于它们之间比公式1、3略高。(3)数值积分计算值【问题】利用数值积分计算，分别用“梯形法”和“Simpson法”精确到10位有效数字，再用“Simpson法”精确到15位有效数字.【解】在本题中，“梯形法”和“Simpson法”的计算原理如下：为利用数值积分计算值，现编写如下两个M文件： function PI=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);function PI

14、=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)梯形法：function m2_19(m) %该函数利用“梯形法”计算指定有效数字位数的，以及最少的划分份数nn0=1;PI=vpa(pi,100);while abs(PI-m2_17

15、(n0)=10(-m+1) n0=2*n0;enda=n0/2;b=n0;n=(a+b)/2; %此处采用了“二分法”的思想，以便于快速找出nwhile b-a=10 %本应为b-a=1，此处是为了避免四舍五入产生的误差导致计算次数过多 if abs(PI-m2_17(n)=10(-m+1) a=n; n=floor(a+b)/2); else b=n; n=floor(a+b)/2); endendwhile abs(PI-m2_17(n)=10(-m+1) n=n+1;enddisp(划分份数至少为n=,num2str(n)vpa(m2_17(n),m)(2)Simpson法： funct

16、ion m2_20(m) %该函数利用“Simpson法”计算指定有效数字位数的，以及最少的划分份数nn=1;PI=vpa(pi,100);while abs(PI-m2_18(n)=10(-m+1) n=n+1;enddisp(划分份数至少为n=,num2str(2*n)vpa(m2_18(n),m)运行结果见第14页图3.1.图3.1【答】由图3.1知，利用“梯形法”和“Simpson法”计算10位有效数字的值所需的最少划分份数分别为12910和20；利用“Simpson法”计算15位有效数字的值所需的最少划分份数为92.此外不难发现，当指定相同有效数字位数时，“Simpson法”所需的划

17、分份数远少于“梯形法”. 换言之，当划分份数一定时，“Simpson法”的计算精度远高于“梯形法”. 综上所述，“Simpson法”的计算效率远高于“梯形法”.（4）利用Monte Carlo法计算值【问题】(1)试用Monte Carlo法计算的5位有效数字； (2)设计方案试用计算机程序模拟Buffon实验.【解】首先将Monte Carlo法和Buffon实验的原理进行一些简要说明.function PI=m2_21(n) %Monte Carlo投点法m=0;for k=1:n if rand2+rand2=1 m=m+1; end;end;PI=vpa(4*m/n,10); 由于取n

18、=10000运行该程序1次和取n=1000运行该程序10次再求平均值的本质没有区别，故每次进行投点实验时，仅运行该程序一次. 运行结果如下：图4.1从图4.1中可以看出，随着投点次数不断增加，所得的值与真实值的误差总体上不断减小. 本次实验中，当投点次数为10亿，可得的5位有效数字. 然而，投点10万次，即可求出的4位有效数字；而投点1000万次，仅可求出的3位有效数字. 这说明增加投点次数，不能说计算的精度必定提高.图4.2图4.3为利用计算机程序进行随机模拟，现编写M文件如下：function m2_29(l,d,m) %该函数利用Buffon投针法计算的近似值%其中l为针长，d为平行线间

19、距，m为试验次数if l/d1 disp(错误！请重新输入l，d值。)else n=0; for i=1:m theta=pi*rand(1); h=-d/2+d*rand(1); if abs(h)=10(-m+1) n=n+1; if mod(n,2)=1 a(n)=(n+1)/n; else a(n)=n/(n+1); end S=S*a(n);enddisp(为达到,num2str(m),位有效数字，所需的项数n=,num2str(n)vpa(2*S,m)运行结果如下：图5.1【答】图5.1表明，利用“Wallis公式”计算值，每增加一位有效数字，所需的项数n就要提高一个数量级. 因此

20、，“Wallis公式”计算的效率较低，大致和反正切函数公式的幂级数展开式的计算效率处于同一数量级.（6）e的计算【问题】试利用多种方法计算自然对数的底数e的近似值.【解】以下四种方法均可用于计算e值.(1)利用公式1：.function m2_24(m) %利用(1+1/n)n计算e的近似值E，m为有效数字位数E=exp(1);n=1;while abs(E-(1+1/n)n)=10(-m+1) n=n+1;enddisp(为达到,num2str(m),位有效数字，所需的n值为,num2str(n)e=vpa(1+1/n)n,m) 证：这是e的定义式，故证明从略. 现编写M文件如下：输入以下命

21、令：图6.1.1运行结果如下： 图6.1.2(2)利用公式2：. 现编写M文件如下：function m2_27(m) %利用(1+1/n)(n+1)计算e的近似值E，m为有效数字位数E=exp(1);n=1;while abs(E-(1+1/n)(n+1)=10(-m+1) n=n+1;enddisp(为达到,num2str(m),位有效数字，所需的n值为,num2str(n)e=vpa(1+1/n)(n+1),m)输入以下命令： 上图：图6.2.1运行结果如下： 下图：图6.2.2(3)利用公式3： .function m2_25(m) %利用Taylor级数计算e的近似值,m为有效数字位

22、数E=vpa(exp(1),100);e=0;n=0;while abs(E-e)=10(-m+1) e=e+1/factorial(n); n=n+1;enddisp(为达到,num2str(m),位有效数字，所需Taylor级数项数n=,num2str(n)e=vpa(e,m) 现编写M文件如下：输入以下命令： 上图：图6.3.1运行结果如下： 下图：图6.3.2(4)利用公式4： 现编写M文件如下：function m2_26(m,k) %利用(1+1/n+k/n2)n计算e的近似值E，m为有效数字位数，k为给定的正常数E=exp(1);n=1;while abs(E-(1+1/n+k/

23、n2)n)=10(-m+1) n=n+1;enddisp(为达到,num2str(m),位有效数字，所需的n值为,num2str(n)e=vpa(1+1/n+k/n2)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的计算效率更高.观察图

24、6.4.3可以发现：对于公式4，当m较小时，参数k越大，计算效率越低（对于较大的m值，由于计算机性能等原因，难以进行实验）.3、实验心得通过本次实验，我对于的几类常见计算方法及其计算效率有了初步的了解，对于MATLAB等编程软件的使用也不再陌生，同时也复习回顾了所学的微积分及概率等知识，感觉收获颇丰.然而，在本次实验中，我也遭遇了一个问题，即所谓的“MATLAB自动补全功能”（参见第6页中部）. 为此，我做了两个小实验进行检验.小实验1:利用的级数展开计算，取n=20. 编写M文件如下：function S=m2_14(n)S=0;for i=1:n if mod(i,2)=0 S=S-1/(

25、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); endendS=vpa(4*S,30); %观察30位有效数注：上述M文件与数学实验（第二版）第22页下部的程序完全相同.运行结果如下：奇怪的是，相同的程序却得出了不同的结果. 更令人意外的是，ans与的前30位竟完全吻合，我认为这可能是MATLAB升级后新增了“自动补全功能”.小实验2：function PI=m2_15(n) %利用Ramanujan公式计算的近似值S=0;for i=0:n S=S+factorial(4*i)*(1103+26390*i)/(factorial(i)4*396(4*i);endS=2*sqrt(2)/9801*S;PI=vpa(1/S,100);利用拉马努金公式计算值，编写M文件：运行结果如下：令人震惊的是，借助拉马努金公式，仅用n=0和n=1两项即得到了的100位有效数字，这绝不可能. 因此我猜测，这应该是MATLAB的“自动补全功能”.不论所谓的“自动补全功能”是否存在，这一现象直接导致在第二部分实验中，当要求的有效数字位数m较大时，利用MATLAB无法得出正确的所需项数n.