电子科技大学《数值计算方法》
实
验
报
告
s1+=(*f)(temp1);/*公式(22)*/
s2+=(*f)(temp2);/*公式(22)*/ }
s=h/6.f*((*f)(a)+4*s1+2*s2+(*f)(b));/*复化辛卜生公式*/
printf("S(%d)=%11.10f\n",n,s);
fprintf(fp,"S(%d)=%11.10f\n",n,s);
deta=s-s0; /*前后两值的间距*/
if(fabs(deta)>eps)/*判断前后所求值的间距是否已满足精度*/
{ n++;/*区间等分数加1*/
s1=0;
s2=0;
s0=s; }
else
break;
}
○3数据处理与分析
当a=0,b=1,eps=1e-6;时
S(1)=0.7833333333，S(2)=0.7853921569
S(3)=0.7853979452，S(4)=0.7853981256
从计算结果看n=4时满足精度要求，即n=4时达到6为有效数字。
当a=0,b=1,eps=5e-10时
S(1)=0.7833333333，S(2)=0.7853921569，S(3)=0.7853979452
S(4)=0.7853981256，S(5)=0.7853981535，S(6)=0.7853981601
S(7)=0.7853981621，S(8)=0.7853981628，S(9)=0.7853981631
从计算结果看n=9时满足精度要求。
在n=7时有9位有效数字，而同时可知n=4时也已达到7位有效数字。
MATLAB绘图
temp2=a+h/2.;
c2=(*f)(temp2);/*公式(25)在k=0下的值*/
temp3=a+3*h/4.;
c3=(*f)(temp3);/*公式(26)在k=0下的值*/
for(k=1;k!=n;k++)
temp1=a+k*h+h/4.;
{
temp2=a+k*h+h/2;
temp3=a+k*h+3*h/4;
temp4=a+k*h;
c1+=(*f)(temp1);/*公式(24)*/
c2+=(*f)(temp2);/*公式(25)*/
c3+=(*f)(temp3);/*公式(26)*/
c4+=(*f)(temp4);/*公式(27)*/ }
c=h/90.f*(7*(*f)(a)+32*c1+12*c2+32*c3+14*c4+7*(*f)(b));/*复化复化柯特斯公式*/ printf("C(%d)=%11.10f\n",n,c);
fprintf(fp,"C(%d)=%11.10f\n",n,c);
deta=c-c0;/*前后两值的间距*/
if(fabs(deta)>eps)/*判断前后所求值的间距是否已满足精度*/
{ n++;/*区间等分数加1*/
c1=0;
c2=0;
c3=0;
c4=0;
c0=c; }
else
break;
}
○3数据处理与分析
当a=0,b=1,eps=1e-6;时
t1=t2;
h/=2.0;/*把区间n等分，每个区间的长度*/
}
else
break;
}
○3数据处理与分析
当a=0,b=1,eps=1e-6;时
T(2)=0.7750000000，T(4)=0.7827941176，T(8)=0.7847471236
T(16)=0.7852354030，T(32)=0.7853574733，T(64)=0.7853879909
T(128)=0.7853956203，T(256)=0.7853975276，T(512)=0.7853980045
从计算结果看n=512时满足精度要求，即n=512时达到6位有效数字。
当a=0,b=1,eps=5e-10时
T(2)=0.7750000000，T(4)=0.7827941176，T(8)=0.7847471236
T(16)=0.7852354030，T(32)=0.7853574733，T(64)=0.7853879909
T(128)=0.7853956203，T(256)=0.7853975276，T(512)=0.7853980045
T(1024)=0.7853981237，T(2048)=0.7853981535，T(4096)=0.7853981609
T(8192)=0.7853981628，T(16384)=0.7853981632
从计算结果看n=16384，满足精度要求，在n=4096时达到8位有效数字。
f 12121212,,,,/2,*2r r c c s s t t h n ======
○
2主要程序代码 while(1)
{
temp=0;
for(k=0;k!=n;k++)
{
x=a+k*h+h/2.;/*公式(31)*/
temp+=(*f)(x);/*公式(32)*/
}
t2=(t1+temp*h)/2.f;/*公式(33)*/
s2=t2+(t2-t1)/3;/*公式(34)*/
/*核心算法中c 步*/
if(n>=2)
{
c2=s2+(s2-s1)/15; }
/*核心算法中d 步*/
if(n>=4)
{
r2=c2+(c2-c1)/63;
/*核心算法中e 步*/
printf("R(%d)=%11.10f\n",n/4,r2);
fprintf(fp,"R(%d)=%11.10f\n",n/4,r2);
deta1=r2-r1;
if(fabs(deta1)<eps)break;
}
r1=r2;c1=c2;s1=s2;t1=t2;h/=2.0;n*=2;/*核心算法中f 步*/
}
○3数据处理与分析
当a=0,b=1,eps=1e-6;时
R(1)=0.7853964459，R(2)=0.7853981596，R(4)=0.7853981634
从计算结果看n=4时满足精度要求，而n=2时已有6位有效数字
当a=0,b=1,eps=5e-10时
R(1)=0.7853964459，R(2)=0.7853981596
R(4)=0.7853981634，R(8)=0.7853981634
从计算结果看n=8时满足精度要求，在n=4时已有至少9位有效数字。
4.实验结论：
○1由图7与图8及计算结果在n=45和n=505下都只达到6为有效数字知复化梯形求积公式其逼近真值方式是来回振荡的。
同时在满足精度要求算得的值不一定是有效的。
但从图6知各个计算所得数据位数逐渐稳定，而从数据看知当计算结果某位数与真值相同该位固定知其一般无效位是后几位。
○2从复化求积公式计算结果中对同一积分函数求解在满足精度要求下看，知其收敛最快的是复化柯特斯公式，其次是复化辛卜生公式，最后是复化梯形求积公式。
○3从自适应梯形公式计算结果看，其并没有出现来回振荡现象，且在初次达6位有效数字时以复化辛卜生公式和复化柯特斯公式具有相同的有效数字，再从图10逼近规律看知各个计算所得数据位数逐渐稳定且从各计算结果看知计算某位数与真值相同该位就固定，可见其计算通常不会出现振荡情况。
○4从对同积分函数所用的各方法看知龙贝格算法收敛最快。
报告评分：
指导教师签字：。