《数值计算方法》实验-(6)

电子科技大学《数值计算方法》实验报告学生姓名：陈文龙学号：2704102014教师：赖生建一、实验室名称：211大楼206室二、实验项目名称：《数值计算方法》实验6——数值积分三、实验原理：1、复化求积公式复化求积公式是为提高计算的积分精度，把积分区间分为若干个小区间，将所求积分()If写成这些小区间上的积分之和，然后对每个小区间上的积分应用柯特斯公式，并把每个小区间上的结果累加即所得为复化柯特斯公式1)、复化梯形公式复化梯形公式为110()[()()]2nnkkkhTffxfx−+==+∑或1011()[()2()()]2nnknkTfhfxfxfx−==++∑，计算过程为：○1令1()/,0;hbant=−=           ○2对1,2,,1kn=−计算  11()ttfakh=++ (1)          ○3 1(()2())2htfatfb=++ (2)           图1复化梯形公式求解积分流程图从流程图看第3，4，5，6，7为其主要步骤，是算法的核心部分，其中第5步为复化梯形公式。       2)、复化辛卜生公式复化辛卜生公式为110[()4()()]62nnkkkkhhSfxfxfx−+==+++∑，计算过程为：○1令12()/,(/2),0;hbansfahs=−=+=输入a,b,eps0tt−ε1(()2())2htfatfb=++bahn−=1010nnttt=+==10,1tn==11()1,2,3,1ttfakhkn=++=−输出tYN结束开始○2对1,2,,1kn=−计算 1122(/2),()ssfakhhssfakh=+++=++(3)○312(()42())6hsfassfb=+++(4)图2复化辛卜生公式求解积分流程图从流程图看第3，4，5，6，7为其主要步骤，是算法的核心部分，其中第5步为复化辛卜生公式。输入a,b,eps0ss−ε12(()42())6hsfassfb=+++bahn−=12010nnssss=+===12(/2),0,1sfahsn=+==1122(/2)()1,2,3,1ssfakhhssfakhkn=+++=++=−输出sYN开始结束3)、复化柯特斯公式复化柯特斯公式为113424110()[7()32()12()32()7()]90nnkkkkkkhCffxfxfxfxfx−++++==++++∑，计算过程为：○1令1234()/,(/4),(/2),(3/4),0;hbancfahcfahcfahc=−=+=+=+=○2对1,2,,1kn=−计算 11(/4)ccfakhh=+++(5)22(/2)ccfakhh=+++(6)33(3/4)ccfakhh=+++(7)44()ccfakh=++(8)○31234(7()321232147())90hcfaccccfb=+++++(9)               图3复化柯特斯公式求解积分流程图从流程图看第2步是初始化即第4步k=0的情况。而第3，4，5，6，7是其主要步骤，是算法的核心部分，其中第5步为复化柯特斯公式。输入a,b,eps012341,,000nncccccc=+=====0cc−ε1234(7()321232147())90hcfaccccfb=+++++bahn−=1234(/4),(/2)(3/4),0,1cfahcfahcfahcn=+=+=+==11223344(/4)(/2)(3/4)()1,2,3,1ccfakhhccfakhhccfakhhccfakhkn=+++=+++=+++=++=−输出cYN开始结束2、自适应梯形公式自适应梯形公式是由复化梯形公式演化而来的，由于复化求积公式中的截断误差的表达式可知加密节点，可以提高求积公式的精度，但在使用求积公式之前必须给出合适的步长，这是一个难题。步长取得太大，满足不了精度；步长取得太小，增加了不必要的运算。因此通常采用把区间逐次二分，反复利用求积公式进行计算直至所求得的前后二次积分值满足精度为止。自适应梯形求积公式就是利用复化梯形公式把区间逐次二分直到满足精度为止。具体如下：变步长梯形算法依据公式121201()()22nnnnnkkhbaTTfxhn−+=−=+=∑。计算时可按如下步骤：○1输入精度1,1,,[()()];2banhbaTfafb−ε==−=+○20;s=○3对0,1,2,,1kn=−计算  (),*/2ssfxxakhh=+=++ (10)○4211(*),2;2TThsnn=+= ○5如果21,TT−ε则结束，否则执行○6；○612/2,,hhTT==转○2。图4自适应梯形公式求解积分流程图从流程图看第2步是初始化即第复化梯形公式n=1的情况。而第3，4，5，6，7是其主要步骤，是算法的核心部分，其中第5步为复化梯形公式中n=2n情况。3、龙贝格求积公式龙贝格求积公式继承了复化梯形梯形公式算法简单的优点，形成的新算法。其主要思想12/2hhtt==输入a,b,eps21tt−ε211(*),22tthsnn=+=1,1(()())/2hbanthfafb=−==+(),/20,1,2,3,1ssfxxakhhkn=+=++=−输出t2YN0s=开始结束是用事后误差估计对近似值进行修正，以提高其精度，使用时只需要在对区间逐次二分的过程中，对梯形值进行加权平均，逐步生成精度更高的积分值。其有以下变换：110()[()()]2nnkkkhTffxfx−+==+∑(11)121201()()()22nnnkkhTfTffx−+==+∑(12)221()3nnnnSTTT=+−(13)221()15nnnnCSSS=+−(14)221()63nnnnRCCC=+−(15)算法如下：○1输入a,b,eps,令1,1,(()())/2hbanthfafb=−==+则○20temp=○3对0,1,2,,1kn=−计算 /2xakhh=++(16)()temptempfx=+(17)○4211(*)2tthtemp=+(18)22211()3sttt=+−(19)○5如果2n≥则22211()15csss=+−否则转○8○6如果n≥4则22211()63rccc=+−否则转○8○7输出r2如果21rreps−则结束否则转○8○812121212,,,,/2,*2rrccsstthn======图5龙贝格算法流程图从流程图看而第3，4，5，6，7，9，10，11，12是其主要步骤，是算法的核心部分。12121212,,,/2,*2rrccsstthn======输入a,b,epstemp=01,1,(()())/2hbanthfafb=−==+/2,()0,1,2,3,1xakhhtemptempfxkn=++=+=−211(*)2tthtemp=+2221()/3sttt=+−n=22221()/15csss=+−n=42221()/63rccc=+−输出r221rreps−YYYNNN开始结束四、实验目的：1、通过与实际计算体会各种方法的精确度。2、会编写用龙贝格算法求定积分的程序。五、实验内容： 1、编写用复求积公式求解积分的程序,计算积分12011Ix=+∫，观察n取多少为时达到有效精度，并对三中求积公式进行比较。2、编写自适应梯形公式求解积分的程序，计算积分12011Ix=+∫，观察n取多少为时达到有效精度。 3、编写龙贝格方法求积分程序，计算积分12011Ix=+∫六、实验结果及讨论：1.复化求积公式1)复化梯形求积公式○1核心算法a对1,2,,1kn=−计算  11()ttfakh=++ (20)b1(()2())2htfatfb=++(21)○2主要程序代码while(1){h=(b-a)/n;/*把区间n等分，每个区间的长度*/for(k=1;k!=n;k++){temp1=a+k*h;/*公式(20)*/t1+=(*f)(temp1);/*公式(21)*/}t=h/2.f*((*f)(a)+2*t1+(*f)(b));/*复化梯形求积公式*/printf(T(%d)=%11.10f\n,n,t);fprintf(fp,T(%d)=%11.10f\n,n,t);deta=t-t0;/*前后两值的间距*/if(fabs(deta)eps)/*判断前后所求值的间距是否已满足精度*/{n++;/*区间等分数加1*/t1=0;t0=t;}elsebreak;}○3数据处理与分析当a=0,b=1,eps=1e-6时T(45)=0.7853775873从结果知n=45时计算满足精度要求。用MATLAB绘图05101520253035404500.0050.010.0150.020.0250.03等分区间数nT(n)-T(n-1)复化梯形公式逼近精度曲线图图6复化梯形求积公式求解积分逼近精度曲线图从图可知复化梯形公式逼近速度是逐渐降低的，即各个计算所得数据位数逐渐稳定。当a=0,b=1,eps=5e-10时T(504)=0.7853979994，T(505)=0.7853980000，T(506)=0.7853980007T(507)=0.7853980013，T(508)=0.7853980019，T(509)=0.7853980026T(510)=0.7853980032，T(511)=0.7853980038，T(512)=0.7853980045达到精度要求时n=551,T(551)=0.7853980262。从计算结果中可看出当n=505时有用6位有效数字前次计算是n=45时有6位有效数字，可见复化梯形公式逼近精确值过程是来回振荡的。MATLAB绘图010020030040050060000.0050.010.0150.020.0250.03等分区间数nT(n)-T(n-1)复化梯形公式逼近精度曲线图图7复化梯形求积公式求解积分逼近精度曲线图取400至551中151点重新作图得40045050055045678910111213x10-10等分区间数nT(n)-T(n-1)复化梯形公式逼近精度曲线图图8复化梯形求积公式求解积分部分点逼近精度曲线图从图可知曲线往复振荡说明复化梯形公式在达到精度要求时是在真值附近来回振荡。2)复化辛卜生求积公式○1核心算法a对1,2,,1kn=−计算  1122(/2),()ssfakhhssfakh=+++=++ (22)b 12(()42())6hsfassfb=+++ (23)  ○2主要程序代码while(1){h=(b-a)/n;/*把区间n等分，每个区间的长度*/temp1=a+h/2.;s1=(*f)(temp1);for(k=1;k!=n;k++){temp1=a+k*h+h/2.;temp2=a+k*h;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;}elsebreak;}○3数据处理与分析当a=0,b=1,eps=1e-6;时S(1)=0.7833333333，S(2)=0.7853921569S(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.7853979452S(4)=0.7853981256，S(5)=0.7853981535，S(6)=0.7853981601S