romberg.m

% 龙贝格求积公式
% 注释的地方需要修改
clear

a=4; b=8;                % a为积分下界，b为积分上界
ep=1e-8;                 % 精度要求，1e-n表示10的-n次幂
syms x
f=(4-sin(x)^2)^0.5;      % 函数表达式

h=b-a; 
n=1; 
i=1; 
T(i,1)=h/2*(subs(f,x,a)+subs(f,x,b));
while 1
    h=h/2;
    xi=a+(2*(1:n)-1)*h;
    T(i+1,1)=1/2*T(i,1)+h*sum(subs(f,x,xi));
    for j=1:i
        T(i+1,j+1)=(4^j*T(i+1,j)-T(i,j))/(4^j-1);
    end
    n=2*n;
    if abs(T(i+1,i+1)-T(i,i))<=ep
        break
    end
    i=i+1;
end

vpa(T(i+1,i+1),10)       % 结果有效数字位数