romberg.asv

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

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