function s=rombg(a, b, TOL)
% a 为积分下限
% b 为积分上限
% TOL 为误差容限
% s 为最后Romberg积分的值
% rombg_f为积分函数
n = 1;
h = b - a;
delt = 1; % 误差初值
x = a;
k = 0;
R = zeros(4, 4);
R(1, 1) = h / 2 * (rombg_f(a) + rombg_f(b));
while delt > TOL
% 如果两次计算的误差大于给定误差则进入循环
k = k + 1;
h = h / 2;
s = 0;
for j = 1:n
x = a + h * (2 * j - 1);
s = s + rombg_f(x);
end
R(k+1, 1) = R(k, 1) / 2 + h * s;
n = 2 * n;
for i = 1:k
R(k+1, i+1) = ((4^i) * R(k+1, i) - R(k, i)) / (4^i - 1);
end
delt = abs(R(k+1, k) - R(k+1, k+1)); % 计算前后两次的值
end
s = R(k+1, k+1);
function f=rombg_f(x)
f = x / (4 + x^2);
end
% s = rombg(f, 0, 1.5, 1.e-6)
s = rombg(0, 1.5, 1.e-6)
%=> s = 0.3241