龙贝格算法在MATLAB的实现

一、实验内容

用龙贝格算法计算积分

二、程序清单与运行结果 

M文件代码如下：

function I=romberg(fun,a,b,e)
% 使用龙贝格(Romberg数值求解公式)
% 例如：
% I=romberg(@(x)x^(3/2),0,1,0.000001)
%
% T =
%
%     1 至 3 列
%  
%     0.500000000000000                   0                   0
%     0.426776695296637   0.402368927062182                   0
%     0.407018110857901   0.400431916044989   0.400302781977176
%     0.401812464799974   0.400077249447332   0.400053605007488
%     0.400463401302048   0.400013713469406   0.400009477737544
%     0.400117671209778   0.400002427845688   0.400001675470774
%     0.400029739862529   0.400000429413446   0.400000296184629
%     0.400007491908991   0.400000075924478   0.400000052358547
%  
%     4 至 6 列
%  
%                     0                   0                   0
%                     0                   0                   0
%                     0                   0                   0
%     0.400049649817493                   0                   0
%     0.400008777304688   0.400008617020324                   0
%     0.400001551625270   0.400001523289272   0.400001516355028
%     0.400000274291198   0.400000269282045   0.400000268056232
%     0.400000048488292   0.400000047602790   0.400000047386095
%
%     7 至 8 列
%  
%                     0                   0
%                     0                   0
%                     0                   0
%                     0                   0
%                     0                   0
%                     0                   0
%     0.400000267751397                   0
%     0.400000047332207   0.400000047318753
%
%
% I =
%
%     0.400000047318753
 
% 判断输入参数是否足够
if nargin~=4
    error('请输入需要求积分的f、上界和下界以及误差e')
end
 
k=0; % 迭代次数
n=1; % 区间划分个数
h=b-a; %上下界间距
T(1,1)=h/2*(fun(a)+fun(b));
d=b-a; %误差间距
while e<=d
    k=k+1;
    h=h/2;
    sum=0;
    % 计算第一列T
    for i=1:n
        sum=sum+fun(a+(2*i-1)*h);
    end
    T(k+1,1)=T(k)/2+h*sum;
    % 迭代
    for j=1:k
        T(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j))/(4^j-1);
    end
    n=n*2;
    d=abs(T(k+1,k+1)-T(k,k));
end
T
I=T(k+1,k+1);

运行结果如下：

I=romberg(@(x)x^(3/2),0,1,0.000001)

T =

  1 至 3 列

   0.500000000000000                   0                   0
   0.426776695296637   0.402368927062182                   0
   0.407018110857901   0.400431916044989   0.400302781977176
   0.401812464799974   0.400077249447332   0.400053605007488
   0.400463401302048   0.400013713469406   0.400009477737544
   0.400117671209778   0.400002427845688   0.400001675470774
   0.400029739862529   0.400000429413446   0.400000296184629
   0.400007491908991   0.400000075924478   0.400000052358547

  4 至 6 列

                   0                   0                   0
                   0                   0                   0
                   0                   0                   0
   0.400049649817493                   0                   0
   0.400008777304688   0.400008617020324                   0
   0.400001551625270   0.400001523289272   0.400001516355028
   0.400000274291198   0.400000269282045   0.400000268056232
   0.400000048488292   0.400000047602790   0.400000047386095

  7 至 8 列

                   0                   0
                   0                   0
                   0                   0
                   0                   0
                   0                   0
                   0                   0
   0.400000267751397                   0
   0.400000047332207   0.400000047318753


I =

   0.400000047318753

三、实验总结

在本次实验中，我使用了函数句柄来求值的方法，来避免直接代入求值而因舍入误差导致结果不理想。