romber.m

来自「余向阳博士Matlab教程的源代码。 包括：1.非线性方程的解法；2.线性方程」· M 代码 · 共 30 行

function[R,quad,err,h]=romber(f1441,a,b,n,tol)
%  f是被积函数。
%  a, b分别为是积分的上下限。
%  n＋1是T数表的列数。
%  tol是允许误差。
%  R是T数表。
%  quad是所求积分值。
M=1;
h=b-a;
err=1;
J=0;
R=zeros(4,4);
R(1,1)=h*(feval('f1441',a)+feval('f1441',b))/2;
while((err>tol)&(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('f1441',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)