%
% 梯形公式的加速--外推技巧: Romberg 算法
%
clc;
clear all; 
close all;

f = @(x) sin(x)./x;
a = 0;   b = 1; 
tol = 1e-7;    % 精度要求
IterMax = 5;   % 最大递归层数

T(1,1) = (b-a)/2 * (1 + f(b)); % f(0) = 1
fprintf(' T(1,1)=%.9f\n',T(1,1));

for k = 1 : IterMax    
    h = (b-a)/(2^(k-1));
    T(k+1,1) = T(k,1)/2 + (h/2) * sum(f(a+0.5*h+[0:2^(k-1)-1]*h));
    fprintf(' T(%d,1)=%.9f\t',k+1,T(k+1,1));
    for j = 1 : k
        T(k+1,j+1) = (4^j * T(k+1,j) - T(k,j) )/(4^j-1);
        fprintf(' T(%d,%d)=%.9f\t',k+1,j+1,T(k+1,j+1));
    end
    fprintf('\n')
    if (abs(T(k+1,k+1)-T(k,k))<tol), break, end        
end


