关于Matlab中龙贝格代码的实现

新世纪社畜战士

function [s,k]=Romberg_li(a,b,eps)
if nargin==2
    eps=1.0e-6;
elseif nargin<2
    error
    return
end
t1=10^4;
t2=-10^4;n=2;
t(1,1)=0.5*(b-a)*(f(a)+f(b));area=t(1,1)/2;
while abs(t2-t1)>=eps
    h=(b-a)/2^(n-1);
    for i=1:22^(n-1))
        area=area+h*f(h*i+a);
    end
    t(n,1)=area ;
    for j=2:n
        for i=1:(n-j+1)
            t(i,j)=(4^(j-1)*t(i+1,j-1)-t(i,j-1))/(4^(j-1)-1);
        end
    end
    t1=t(i,n);
    t2=t(i+1,n-1);
    n=n+1;
    area=area/2;
end
s=t1;k=n-1;