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.
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.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
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);
matrix(ii,1) = matrix(ii,1) + g(b);
matrix(ii,1) = matrix(ii,1)*h/2;
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);
disp(‘Number of Layers required to be an Integer’)
matrix = matrix
I = matrix(n,n);