romberg.m

为配套王能超版本数值分析简明教程编写,方便学生教师演示!

function X=romberg(f1,a,b,n)
l=1;
if n>=4
    f=inline(f1);
    h=b-a;
    m=1;
    for i=1:n
        s1=0;
        for k=1:m-1
            s1=f(a+h*k)+s1;
        end
        T(i)=h/2*(f(a)+f(b))+h*s1;
        h=h/2;
        m=m*2;
        if i>=2
            S(i-1)=4/3*T(i-1)-1/3*T(i);
        end
        if i>=3
            C(i-2)=-1/15*S(i-2)+16/15*S(i-1);
        end
        if i>=4
            R(i-3)=-1/63*C(i-3)+64/63*C(i-2);
            fprintf('%d%13.8f\n',l,R(i-3));
            l=l+1;
        end
    end
    X=R(n-3);
else
    disp('n必须为大于3的整数');
end