Romberg Integration in Numerical Analysis

romberg.m

% Romberg Integration in Numerical Analysis
% Author: Farhad Mehryari
% Web: blog.mrfarhad.ir
function r = romberg(f,a,b,n)
h = (b - a) ./ (2.^(0:n-1));
r(1,1) = (b - a) * (f(a) + f(b)) / 2;
for j = 2:n
    subtotal = 0;
    for i = 1:2^(j-2)
        subtotal = subtotal + f(a + (2 * i - 1) * h(j));
    end
    r(j,1) = r(j-1,1) / 2 + h(j) * subtotal;
    for k = 2:j
        r(j,k) = (4^(k-1) * r(j,k-1) - r(j-1,k-1)) / (4^(k-1) - 1);
    end
end;
% Author: Farhad Mehryari
% Usage: نحوه ی استفاده
r=romberg(@sin,0,pi,5)
r =

   0.00000   0.00000   0.00000   0.00000   0.00000
   1.57080   2.09440   0.00000   0.00000   0.00000
   1.89612   2.00456   1.99857   0.00000   0.00000
   1.97423   2.00027   1.99998   2.00001   0.00000
   1.99357   2.00002   2.00000   2.00000   2.00000
% Enjoy :)