数值分析实验报告.doc

数值分析上机实验报告数值分析实验报告实验名称复化辛普森公式求定积分实验时间2009.10姓名班级学号成绩一、实验目的，内容 1理解复化梯形公式、复化Simpson公式、Romberg方法和复化Gauss-Legendre公式计算的概念2掌握Newton-Cotes求积公式的原理，包括了解这些公式的误差及代数精度，参考课本写出复化辛普森算法的程序，在matlab中实现，并用matlab内置的函数计算，进行误差分析。三代码 function y=simpson(f,a,b,yxsz)if narginp) for j=1:2 n=n*j; h=(b-a)/(2*n); s1=0.0;s2=0.0; for i=1:n-1 x=a+2*i*h; s1=s1+eval(f); end for i=1:n x=a+(2*i-1)*h; s2=s2+eval(f); end sn(j)=h*(subs(f,a)+4*s2+2*s1+subs(f,b)/3; end fprintf(n=%3d Sn=%.10fn,n,sn(1); eps=abs(sn(1)-sn(2);endtocny=sn(2);程序说明： 被积函数f可在调用simpson时作为一个参数输入，是字符串型的函数，程序中用eval计算其在某特定点的值(用subs也可).最外侧while每次循环时先计算区间等分数为n的积分值sn(1)，然后将n加倍，再计算一次sn(2)，若两次计算结果之差小于指定精度p，while的循环条件将不满足，停止计算，此时的sn(2)作为积分近似值。一、 实验结果 y=simpson(1/(1+x2),0,1,6)n= 2 Sn=0.7833333333n= 4 Sn=0.7853921569n= 8 Sn=0.7853981256Elapsed time is 0.453000 seconds.n = 8y = 0.785398162806206四.数值结果 y=simpson(1/(1+x2),0,1,6)n= 2 Sn=0.7833333333n= 4 Sn=0.7853921569n= 8 Sn=0.7853981256Elapsed time is 0.453000 seconds.n = 8y = 0.785398162806206五计算结果的分析在fun.m中定义函数：function y = fun(x) y = 1./(x.2+1);在命令窗口中调用：quad(fun,0,1)ans = 0.78539814924326quad命令也是用复化辛普森方法，默认精度是1e-6。上述积分的精确值是pi/4= 0.785398163397448六. 计算中出现的问题，解决方法及体会区间等分数从1开始逐次加倍，若两次相差不超过误差限，则停止计算。在计算的过程中每次都输出区间等分数和对应的积分值，这样做是为了观察收敛的速度和计算过程，但明显减慢了计算速度。从理论上分析可知，如果被积函数不满足光滑性的要求，那么辛普森算法的效率也会大打折扣。如计算：y=simpson(sqrt(x),0,1,6)n= 2 Sn=0.6380711875n= 4 Sn=0.6565262648n= 8 Sn=0.6630792801n= 16 Sn=0.6653981886n= 32 Sn=0.6662181827n= 64 Sn=0.6665081031n=128 Sn=0.6666106059n=256 Sn=0.6666468462n=512 Sn=0.6666596591n=1024 Sn=0.6666641891n=2048 Sn=0.6666657907Elapsed time is 3.406000 seconds.n = 2048y = 0.666666356971916这个积分一眼就能看出是2/3,但辛普森方法却用了3.4秒将区间划分了2048等分才达到6位有效数字。用quad命令： quad(fun,0,1)ans = 0.66665956927202Elapsed time is 0.000000 seconds.其中fum.m中定义了：function y = fun(x