romber.m

用于求积分

%龙贝格求积公式
function [R,quad,err,h]=romber(f,a,b,n,delta)
%f被积函数
%a，b是上下限
%n+1是T数表的列数
%delta是允许误差
%R是T的数表
%quad是所求的积分
M=1;
h=b-a;
err=1;
J=0;
R=zeros(4,4);
R(1,1)=h*(feval(f,a)+feval(f,b))/2;
while((err>delta)&(J<n))|(J<4)
    J=J+1;h=h/2;s=0;
    for p=1:M
        x=a+h*(2*p-1);
        s=s+feval(f,x);
    end
    R(J+1,1)=R(J,1)/2+h*s;
    M=2*M;
    for K=1:J
        R(J+1,K+1)=R(J+1,K)+(R(J+1,K)-R(J,K))/(4^K-1);
    end
    err=abs(R(J,J)-R(J+1,K+1));
end
quad=R(J+1,J+1);