Topic 13.3: Romberg Integration (Matlab)

Introduction Notes Theory HOWTO Examples Engineering Error Questions Matlab Maple

Approximate the integral of f(x) = cos(x) until ε_step < 1e-10 or N = 10.

Please note, because we cannot begin at R(0, 0), each index will be off by one. This is not optimal code, either.

format long
eps_step = 1e-10;
N = 10;
R = zeros( N + 1, N + 1 );
a = 0;
b = 10;
h = b - a;
R(0 + 1, 0 + 1) = 0.5*(cos(a) + cos(b))*h;

for i = 1:N
    h = h/2;
    % This calculates the trapezoidal rule with intervals of width h
    R( i + 1, 1 ) = 0.5*(cos(a) + 2*sum( cos( (a + h):h:(b - h) ) ) + cos(b))*h;

    for j = 1:i
        R(i + 1, j + 1) = (4^j*R(i + 1, j) - R(i, j))/(4^j - 1);
    end

    if abs( R(i + 1, i + 1) - R(i, i) ) < eps_step
       break;
    elseif i == N + 1
       error( 'Romberg integration did not converge' );
    end
end

R( i + 1, i + 1 )

The correct answer to 20 decimal digits is sin(10) = -0.54402111088936981340.