lomberg.m

来自「数值分析中的插值计算」· M 代码 · 共 36 行

function romberg(a,b,n,tol)
%Romberg积分
%Input   
%         a,b是积分下限与上限
%         n 是设定的最大的区间等分数
%         tol 限差
%Output   R Romberg表
%         quad 求积值
%         err  误差估计
syms x;
f=input('f=');% f is the function
M=1;
h=b-a;
err=1;
J=0;
R=zeros(4,4);
R(1,1)=h*(subs(f,a)+subs(f,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+subs(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);
R
quad
err