利用Matlab 实现Romberg 数值积分算法
一、内容摘要
针对于某些多项式积分，利用Newton —Leibniz 积分公式求解时有困难，可以采用数值积分的方法，求解指定精度的近似解，本文利用Matlab 中的.m 文件编写了复化梯形公式与Romberg 的数值积分算法的程序，求解多项式的数值积分，比较两者的收敛速度。
二、数值积分公式
1.复化梯形公式求解数值积分的基础是将区间一等分时的Newton —Cotes
求积公式:
I =(x)[f(a)f(b)]2
b
a b a
f dx -≈
+⎰ 其几何意义是，利用区间端点的函数值、与端点构成的梯形面积来近似(x)f 在区间[a,b]上的积分值，截断误差为：
3"
(b a)()12
f η-- (a,b)η∈ 具有一次的代数精度，很明显，这样的近似求解精度很难满足计算的要求，因而，可以采用将积分区间不停地对分，当区间足够小的时候，利用梯形公式求解每一个小区间的积分近似值，然后将所有的区间加起来，作为被求函数的积分，可以根据计算精度的要求，划分对分的区间个数，得到复化梯形公式：
I =1
1
(b a)(b a)
(x)dx [f(a)f(b)2(a )]2n b
a
k k f f n n -=--≈+++∑⎰
其截断误差为：
2"
(b a)h ()12
R f η--=
(a,b)η∈ 2.Romberg 数值积分算法
使用复化的梯形公式计算的数值积分，其收敛速度比减慢，为此，采用Romberg 数值积分。
其思想主要是，根据I 的近似值2n T 加上I 与2n T 的近似误差，作为新的I 的近视，反复迭代，求出满足计算精度的近似解。
用2n T 近似I 所产生的误差可用下式进行估算：
12221
()3
n n n I T T T -∆=-=-
新的I 的近似值：
122
n n j T T -=∆+ j =（0 1 2 ….） Romberg 数值积分算法计算顺序
i=0 (1) 002T
i=1 (2) 102T
(3) 012T
i=2 (4) 202T (5) 112T
(6) 022T
i=3
(7)
302T (8) 212T (9) 122T
(10) 032T i=4 (11)
402T
(12) 312T
(13) 222T
(14) 132T
…
…
…
…
其中，第一列是二阶收敛的，第二列是四阶收敛的，第三列是六阶收敛的，第四列是八阶收敛的，即Romberg 序列。
三、复化梯形法以及Romberg算法程序流程图
图1 复化梯形法程序流程图
图2 Romberg算法程序流程图
四、计算实例
依据上文所述的流程图，编写复化梯形程序以及Romberg 算法程序，并且利用实例验证程序的正确性，示例如下(计算精度)：
1
2
04
1dx x π=+⎰
表2 计算结果
计算精度 0.5×10^-5 0.5×10^-7
0.5×10^-9
复化梯形 算法 时间
0.069826394633 0.216635802304 3.459824945493
近似值 3.141590110458 3.141592613853 3.141592653434 Romberg 算法
时间
0.0456******** 0.0433******** 0.044913907518
近似值 3.141592502458 3.141592651224 3.141592653552
从上表中可以看出，当要求的计算精度不高时，复化梯形算法与Romberg 算法计算时间相差不太大，但是Romberg 算法是要快于复化梯形算法的；当要求的计算精度更高的时候，Romberg 算法是明显快于复化梯形算法。
本文所编写的程序适用于多项式的数值积分，且对于积分区间内，被积函数在每一点必须有定义，在以后的学习中进一步改进。
附录：
1.复化梯形算法程序
function []=sf(a,b,m,M,d)
tic
disp('请输入分子多项式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^30
t=0;
for j=0:1:2^(i-1)-1
v=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 %存储其他项
f=i+1; %记录符合条件的值的下标
if(1/3*(kk(i+1,1)-kk(i,1))<=d)
break;
end
end
time=toc
fprintf('The result is %f\n', kk(f,1))
2.Romberg算法程序
function []=romberg(a,b,m,M,d)
tic
disp('请输入分子多项式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^40
t=0;
for j=0:1:2^(i-1)-1
v=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)))<=d) %判断梯形公式值是否达到要求disp('The result is:')
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=(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('The result is:')
kk(i+1,3)
break;
else
r=(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('The result is:')
kk(i+1,4)
break;
end
end
end
end
end
end
end
end
time=toc。