利用Matlab实现Romberg数值积分算法----系统建模与仿真结课作业

利用Matlab实现Romberg数值积分算法一、内容摘要针对于某些多项式积分，利用Newton—Leibniz积分公式求解时有困难，可以采用数值积分的方法，求解指定精度的近似解，本文利用Matlab中的.m文件编写了复化梯形公式与Romberg的数值积分算法的程序，求解多项式的数值积分，比较两者的收敛速度。二、数值积分公式1.复化梯形公式求解数值积分的基础是将区间一等分时的Newton—Cotes求积公式:I=(x)[f(a)f(b)]2babafdx其几何意义是，利用区间端点的函数值、与端点构成的梯形面积来近似(x)f在区间[a,b]上的积分值，截断误差为：3(ba)()12f(a,b)具有一次的代数精度，很明显，这样的近似求解精度很难满足计算的要求，因而，可以采用将积分区间不停地对分，当区间足够小的时候，利用梯形公式求解每一个小区间的积分近似值，然后将所有的区间加起来，作为被求函数的积分，可以根据计算精度的要求，划分对分的区间个数，得到复化梯形公式：I=11(ba)(ba)(x)dx[f(a)f(b)2(a)]2nbakkffnn其截断误差为：2(ba)h()12Rf(a,b)2.Romberg数值积分算法使用复化的梯形公式计算的数值积分，其收敛速度比减慢，为此，采用Romberg数值积分。其思想主要是，根据I的近似值2nT加上I与2nT的近似误差，作为新的I的近视，反复迭代，求出满足计算精度的近似解。用2nT近似I所产生的误差可用下式进行估算：12221()3nnnITTT新的I的近似值：122nnjTTj=（012….）Romberg数值积分算法计算顺序i=0(1)002Ti=1(2)102T(3)012Ti=2(4)202T(5)112T(6)022Ti=3(7)302T(8)212T(9)122T(10)032Ti=4(11)402T(12)312T(13)222T(14)132T…………其中，第一列是二阶收敛的，第二列是四阶收敛的，第三列是六阶收敛的，第四列是八阶收敛的，即Romberg序列。三、复化梯形法以及Romberg算法程序流程图输入被积函数、积分区间、收敛条件计算Ti0计算Ti+10差值是否满足收敛条件输出结果是否图1复化梯形法程序流程图输入被积函数、积分区间、收敛条件计算并存储Ti0计算并存储Ti+10差值是否满足收敛条件输出结果是构建并存储T1i+1i是否大于等于2否是差值是否满足收敛条件是否否构建并存储T2ii是否大于等于3否是差值是否满足收敛条件是否构建并存储T3ii是否大于等于4否是差值是否满足收敛条件否是图2Romberg算法程序流程图四、计算实例依据上文所述的流程图，编写复化梯形程序以及Romberg算法程序，并且利用实例验证程序的正确性，示例如下(计算精度)：12041dxx表2计算结果计算精度0.5×10^-50.5×10^-70.5×10^-9复化梯形算法时间0.0698263946330.2166358023043.459824945493近似值3.1415901104583.1415926138533.141592653434Romberg算法时间0.0456873297100.0433617263570.044913907518近似值3.1415925024583.1415926512243.141592653552从上表中可以看出，当要求的计算精度不高时，复化梯形算法与Romberg算法计算时间相差不太大，但是Romberg算法是要快于复化梯形算法的；当要求的计算精度更高的时候，Romberg算法是明显快于复化梯形算法。本文所编写的程序适用于多项式的数值积分，且对于积分区间内，被积函数在每一点必须有定义，在以后的学习中进一步改进。附录：1.复化梯形算法程序function[]=sf(a,b,m,M,d)ticdisp('请输入分子多项式a，分母多项式b，积分下限m,积分上限M,以及计算精度d')f=poly2sym(a)/poly2sym(b)%用于给用户显示被积函数的形式%利用梯形公式计算此数值积分disp('利用梯形公式计算数值积分的结果')kk=zeros();%用于存放结果kk(1,1)=1/2*(M-m)/1*(subs(f,'x',m)+subs(f,'x',M))%先存储首项fori=1:1:2^30t=0;forj=0:1:2^(i-1)-1v=m+(2*j+1)*(M-m)/(2^i)vv=polyval(a,v)/polyval(b,v);t=t+(M-m)/(2^i)*vvendy=1/2*kk(i,1)+t%通项公式计算各项值kk(i+1,1)=y%存储其他项f=i+1;%记录符合条件的值的下标if(1/3*(kk(i+1,1)-kk(i,1))=d)break;endendtime=tocfprintf('Theresultis%f\n',kk(f,1))2.Romberg算法程序function[]=romberg(a,b,m,M,d)ticdisp('请输入分子多项式a，分母多项式b，积分下限m,积分上限M,以及计算精度d')f=poly2sym(a)/poly2sym(b)%用于给用户显示被积函数的形式disp('利用梯形公式计算数值积分的结果')kk=zeros();%用于存放结果kk(1,1)=1/2*(M-m)/1*(subs(f,'x',m)+subs(f,'x',M));%先存储首项fori=1:1:2^40t=0;forj=0:1:2^(i-1)-1v=m+(2*j+1)*(M-m)/(2^i);vv=polyval(a,v)/polyval(b,v);t=t+(M-m)/(2^i)*vv;endy=1/2*kk(i,1)+t;%通项公式计算各项值kk(i+1,1)=y;%存储其他项if(abs(1/3*(kk(i+1,1)-kk(i,1)))=d)%判断梯形公式值是否达到要求disp('Theresultis:')kk()kk(i+1,1)%梯形值满足要求，输出结果break;elses=(4*kk(i+1,1)-kk(i,1))/(4-1);%构造simpson各项kk(i+1,2)=s%存储if(i+1=3)if(i+1=3&abs(1/15*(kk(i+1,2)-kk(i,2)))=d)kk()disp('Theresultis:')kk(i+1,2)%simpson值满足要求，输出结果pan1=0;break;elsec=(4^2*kk(i+1,2)-kk(i,2))/(4^2-1);%构造cotes值kk(i+1,3)=c%存储cotes值if(i+1=4)if(i+1=4&abs(1/63*(kk(i+1,3)-kk(i,3)))=d)disp('Theresultis:')kk(i+1,3)break;elser=(4^3*kk(i+1,3)-kk(i,3))/(4^3-1)%构造romberg值kk(i+1,4)=r%存储romberg值if(i+1=5)if(i+1=5&abs(1/127*(kk(i+1,4)-kk(i,4)))=d)disp('Theresultis:')kk(i+1,4)break;endendendendendendendendtime=toc