sf_romberg.m

基于matlab的龙贝格算法程序

function [T]=SF_Romberg(a,b,f,tol,maxit)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%数值积分的Romberg算法
%
%Oct 31，2007
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

k=1;      %k 记录二分次数
h=b-a;          %当前步长
T=zeros(maxit,maxit);

T(1,1)=(b-a)*(f(a)+f(b))/2;         %准备初值
e=1;          %误差
while (e>tol)&(k<maxit)
    
    h1=h/2;     %二分后的新步长
    k=k+1;
    for i=a+h1:h:b-h1
        T(k,1)=T(k,1)+f(i);
    end                             
    T(k,1)=T(k-1,1)/2+T(k,1)*h1;        %求梯形值
    
    for j=2:k
        T(k,j)=(T(k,j-1)*4^(j-1)-T(k-1,j-1))/(4^(j-1)-1);     %求加速值
    end

    h=h1;
    e=abs(T(k,k)-T(k-1,k-1));      %精度控制
end
e