利用Matlab实现Romberg数值积分算法

利用 Matlab 实现 Romberg 数值积分算法一、内容摘要针对于某些多项式积分，利用 NewtonLeibniz 积分公式求解时 有困难，可以采用数值积分的方法，求解指定精度的近似解，本文利 用Mat lab中的.m文件编写了复化梯形公式与Romberg的数值积分算 法的程序，求解多项式的数值积分，比较两者的收敛速度。二、数值积分公式1. 复化梯形公式求解数值积分的基础是将区间一等分时的NewtonCotes 求积公式:f .b aI = J bf (x)dx uf(a) + f(b) a2其几何意义是，利用区间端点的函数值、与端点构成的梯形面积 来近似f(x)在区间a,b上的积分值，截断误差为：(b a)312f (n)耳 e (a,b)具有一次的代数精度，很明显，这样的近似求解精度很难满足计 算的要求，因而，可以采用将积分区间不停地对分，当区间足够小的 时候，利用梯形公式求解每一个小区间的积分近似值，然后将所有的 区间加起来，作为被求函数的积分，可以根据计算精度的要求，划分 对分的区间个数，得到复化梯形公式：I = Jbf(x)dx 导f(a) + f(b) + 2艺f(a+ 4)a2n- n其截断误差为：R = (b a) h2 f (n )X (a，b)12数值积分算法使用复化的梯形公式计算的数值积分，其收敛速度比减慢，为此，采用Romberg数值积分。其思想主要是，根据i的近似值T加上i与T2n2n的近似误差，作为新的I的近视，反复迭代，求出满足计算精度的近 似解。用t近似I所产生的误差可用下式进行估算：2n1A = I T 二一(T T )232n 2n1新的I的近似值：Tj =A + Tj =(0 1 2 .)2n12nRomberg数值积分算法计算顺序i=0(1)T200i=1 (2)T021(3)T120i=2(4)T022(5)T121(6)T220i=3(7)T0(8)T1(9)T2(10)T323222120i=4(11)T024(12)T123(13)T222(14)T321其中，第一列是二阶收敛的，第二列是四阶收敛的，第三列是六阶收敛的，第四列是八阶收敛的，即 Romberg 序列。三、复化梯形法以及Romberg算法程序流程图图 1 复化梯形法程序流程图图2 Romberg算法程序流程图四、计算实例依据上文所述的流程图，编写复化梯形程序以及 Romberg 算法程序，并且利用实例验证程序的正确性，示例如下（计算精度）：表 2 计算结果兀二 j1dx0 1 + X 2计算精度X10八-5X10八-7X10八-9复化梯形时间算法近 似 3.3.3.值Romberg时间算法近 似 3.3.3.值从上表中可以看出，当要求的计算精度不高时，复化梯形算法与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项for i=1：1：2八30t=0;for j=0：1：2八(i1)1v二m+(2*j+1)*(M-m)/(2i)vv=polyval(a,v)/polyval(b,v);t二t+(M-m)/(2八i)*vvendy=1/2*kk(i,1)+t计算各项值%用于存放结,M) %先存储首% 通项公式% 存储其他f=i+1;%记录符合条kk(i+1,1)=y件的值的下标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*(subs(f,x,m)+subs(f,x,M); %先存储首项 for i=1：1：2八40t=0;for j=0：1：2八(i1)1v二m+(2*j+1)*(M-m)/(2八i)；vv=polyval(a,v)/polyval(b,v);t二t+(M-m)/(2八i)*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)=3)if(i+1=3abs(1/15*(kk(i+1,2)-kk(i,2)=4)if(i+1=4&abs(1/63*(kk(i+1,3)-kk(i,3)=5)if(i+1=5 & abs(1/127*( kk(i+1,4)-kk(i,4)=d )disp(The result is:)endendtime=tocendendendendendendbreak;