Romberg Numerical Integration - Matlab Script

Romberg numerical integration uses trapenzoidal intergration to increase the accuracy of calculating an area. It builds up a matrix using the trapenzoidal rule then extends to calculate a more accurate answer.

The script i have made for the Trapezoidal Rule is at the following link. The script for romberg integration has the trpezoidal script in built when calculating the first column.

>>> Trapezoidal Rule of Numerical Integration - Matlab Scripts

The first column is calculated using the trapenzoidal rule for integration where the number of areas the function is broken up into increases by 2 each time:

Row 1: n = 1, Row 2: n = 2, Row 3: n = 4, Row 4: n = 8...16...32...

As for the other columns, they are calculated by using the values of the column to the left.

R(i,j) = ((4^j)-R(i,j-1) - R(i-1,j-1))/((4^j)-1)

R(1,1) = (4*R(1,0) - R(0,0))/(4-1)

           = (4(0.898904) - 0.888511)/3

           = 0.902368

>> romberg('sin(x)/x',1,3,4)

matrix =

   0.888510987494519                   0                                  0                                  0
   0.898904207160100   0.902368613715294                   0                                  0
   0.901644861268860   0.902558412638446   0.902571065899989                   0
   0.902337806742469   0.902568788567005   0.902569480295576   0.902569455127252

ans =

   0.902569455127252

Continueing down, the most accurate answer is the bottom right one.

Romberg Numerical Integration - Matlab Script

function I = romberg(f_str,a,b,n)
%ROMBERG Romberg Rule integration.
% I = ROMBERG(F_STR, A, B, N) returns the Romberg Integration approximation
% for the integral of f(x) from x=A to x=B, to n layers, where
% F_STR is the string representation of f. Also displays the matrix used to
% calculate the inegral.

matrix = zeros(n);
g = inline(f_str);

if (rem(n,1) == 0)
    for ii = 1:n
        h = (b-a)/(2^(ii-1));
        matrix(ii,1) = matrix(ii,1) + g(a);

        for kk = (a+h):h:(b-h)
            matrix(ii,1) = matrix(ii,1) + 2*g(kk);
        end

        matrix(ii,1) = matrix(ii,1) + g(b);
        matrix(ii,1) = matrix(ii,1)*h/2;
    end
for jj = 2:n
    for ii = jj:n
        matrix(ii,jj) = ((4^(jj-1))*matrix(ii,jj-1)-matrix(ii-1,jj-1))/((4^(jj-1))-1);
    end
end

else
    disp('Number of Layers required to be an Integer')
end
matrix = matrix
I = matrix(n,n);

Other Matlab Scripts I have created include:

>>> Simpson's Rule and Trapezoidal Rule of Numerical Integration - Matlab Scripts

>>> Lagrange Method and Newton Divided Difference Method - Matlab Scripts

>>> Newtons Method of finding Roots - Matlab Script

>>> Bisection Method of finding Roots - Matlab Script

>>> Secant Method of finding Roots - Matlab Script