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

1、利用matlab实现romberg数值积分算法-系统建模与仿真结课作业利用matlab实现romberg数值积分算法一、内容摘要针对于某些多项式积分，利用newtonleibniz积分公式求解时有困难，可以采用数值积分的方法，求解指定精度的近似解，本文利用matlab中的.m文件编写了复化梯形公式与romberg的数值积分算法的程序，求解多项式的数值积分，比较两者的收敛速度。二、数值积分公式1.复化梯形公式求解数值积分的基础是将区间一等分时的newtoncotes求积公式:=其几何意义是，利用区间端点的函数值、与端点构成的梯形面积来近似在区间a,b上的积分值，截断误差为： 具有一次的代数精度，

2、很明显，这样的近似求解精度很难满足计算的要求，因而，可以采用将积分区间不停地对分，当区间足够小的时候，利用梯形公式求解每一个小区间的积分近似值，然后将所有的区间加起来，作为被求函数的积分，可以根据计算精度的要求，划分对分的区间个数，得到复化梯形公式：=其截断误差为： 数值积分算法使用复化的梯形公式计算的数值积分，其收敛速度比减慢，为此，采用romberg数值积分。其思想主要是，根据的近似值加上与的近似误差，作为新的的近视，反复迭代，求出满足计算精度的近似解。用近似所产生的误差可用下式进行估算：新的的近似值： =（0 1 2 .）romberg数值积分算法计算顺序i=0(1)i=1(2)(3)i

3、=2(4)(5)(6)i=3(7)(8)(9)(10)i=4(11)(12)(13)(14)其中，第一列是二阶收敛的，第二列是四阶收敛的，第三列是六阶收敛的，第四列是八阶收敛的，即romberg序列。 三、复化梯形法以及romberg算法程序流程图 图1 复化梯形法程序流程图图2 romberg算法程序流程图四、计算实例依据上文所述的流程图，编写复化梯形程序以及romberg算法程序，并且利用实例验证程序的正确性，示例如下(计算精度)：表2 计算结果计算精度×10-5×10-7×10-9复化梯形算法时间近似值3.3.3.romberg算法时间近似值3.3.3.从上

4、表中可以看出，当要求的计算精度不高时，复化梯形算法与romberg算法计算时间相差不太大，但是romberg算法是要快于复化梯形算法的；当要求的计算精度更高的时候，romberg算法是明显快于复化梯形算法。本文所编写的程序适用于多项式的数值积分，且对于积分区间内，被积函数在每一点必须有定义，在以后的学习中进一步改进。附录：1.复化梯形算法程序function =sf(a,b,m,m,d)ticdisp('请输入分子多项式a，分母多项式b，积分下限m,积分上限m,以及计算精度d')f=poly2sym(a)/poly2sym(b) %用于给用户显示被积函数的形式%利用梯形公式计算

5、此数值积分disp('利用梯形公式计算数值积分的结果')kk=zeros(); %用于存放结果kk(1,1)=1/2*(m-m)/1*(subs(f,'x',m)+subs(f,'x',m) %先存储首项for i=1:1:230 t=0; for j=0:1:2(i-1)-1 v=m+(2*j+1)*(m-m)/(2i) vv=polyval(a,v)/polyval(b,v); t=t+(m-m)/(2i)*vv end y=1/2*kk(i,1)+t %通项公式计算各项值 kk(i+1,1)=y %存储其他项 f=i+1; %记录符合条件的

6、值的下标 if(1/3*(kk(i+1,1)-kk(i,1)<=d) break; endendtime=tocfprintf('the result is %fn', kk(f,1)算法程序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

7、*(subs(f,'x',m)+subs(f,'x',m); %先存储首项for i=1:1:240 t=0; for j=0:1:2(i-1)-1 v=m+(2*j+1)*(m-m)/(2i); vv=polyval(a,v)/polyval(b,v); t=t+(m-m)/(2i)*vv; end y=1/2*kk(i,1)+t; %通项公式计算各项值 kk(i+1,1)=y ; %存储其他项 if(abs(1/3*(kk(i+1,1)-kk(i,1)<=d) %判断梯形公式值是否达到要求 disp('the result is:')

8、kk() kk(i+1,1) %梯形值满足要求，输出结果 break; else s=(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('the result is:') kk(i+1,2) %simpson值满足要求，输出结果 pan1=0; break; else c=(42*kk(i+1,2)-kk(i,2)/(42-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('the result is:') kk(i+1,3) break; else r=(43*kk(i+1,3)-kk(i,3)/(43-1)%构造romberg值 kk(i+1,4)=r %存储romberg值 if(i+1>=5) if(i+1&