I wrote a matlab function to
generate trapezoidal integral for specific number of strips. Then took the
output and called a function I wrote to generate the Romberg table.
This below is the result:
>>
nma_test_romberg
R =
1815.42957071731 665.399801008666 514.07794397635 504.693241461852
952.907243435826 523.53556004087 504.839877438641 0
630.878480889609 506.00835760128 0 0
537.225888423362 0 0 0
The analytical integration
result was 504.53599.
hence, now I can calculate
the required percentage errors:
epsilon(t) is the relative
percentage error compared with the true value. I.e. error normalized to true
value:
true - approx
epsilon(t) = abs
----------------- * 100
true
while epsilon(a) is the error
normalized to current approximate value, i.e.
current step approx - previous step approx
epsilon(a)
= abs ------------------------------------------------- * 100
current approximation
using the above definitions:
I_exact - R(8) 504.53599 – 504.69324146
epsilon(t)
= abs --------------- * 100 = -------------------------- *100
I_excat 504.53599
=
0.0311675 %
Notice also that error in Romberg at column 4 is O(h^8) where h is the strip size. For n=8, h=(3-0)/n = 0.375
Hence error is O (0.000391)
Now for epsilon(a), use for
current approximation I(1,k) and for previous approximation I(1,k-1), where k
is the column number. Last column in the above Romberg table is k=4. Hence I
get:
R(8)– R(6) 504.693241461852 - 514.07794397635
epsilon(a)
= abs ----------- * 100 = abs ------------------------------------- * 100
R(8) 504.693241461852
= 1.85948646 %
The following is the MATLAB code I wrote to help solve this problem.
function nma_test_romberg
%function nma_testRomberg
%
% Tests the romberg interation function nma_romberg
f='x*exp(2*x)';
from=0;
to=3;
nStrips=8;
rowNumber=0;
nSoFar=0;
while(1)
rowNumber=rowNumber+1;
nSoFar=2^(rowNumber-1);
if(nSoFar>nStrips)
break;
end
c(rowNumber)=nma_trapezoidal(f,from,to,nSoFar);
end
R=nma_romberg(c);
R
function R=nma_romberg(firstColumn)
%function nma_romberg(firstColumn)
%
% Function to generate the Romberg integration table given the
% first column (which should be the result of the trapozidal
% integration, for number of strips n=1,2,4,8,16,....
%
% INPUT
% firstColumn: a Vector whose elements represent the integral
% for different number of strips. The first element if trapozidal
% integral for n=1, second element in the vector is for
% trapozidal integral for n=2, the third element is for n=4,
% and so forth. Notice that n goes as 2^i where i is the
% position number in the vector. Use nma_trapezoidal to generate
% the integral for specific n values and use the output to feed
% it to this function.
%
% OUTPUT:
% a Matrix R the reppresents the romberg table. The first column of the
% matrix is the same as the input argument. The remaining columns are
% the result of the extrapolation used.
%
% see nma_testRomberg.m for an example of how to use these functions
%
% Author: Nasser Abbasi
%
% change history
% name date change
% nma 5/2/03 started. To solve HW#4, MAE 185. UCI
%
len=length(firstColumn);
R=zeros(len);
R(:,1)=firstColumn';
nEntries=len;
k=1;
while(nEntries ~= 0)
k=k+1;
F= 4^(k-1);
for(j=1:nEntries-1)
R(j,k)= ( F*R(j+1,k-1) - R(j,k-1) )/ (F-1);
end
nEntries=nEntries-1;
end
function I=nma_trapezoidal(func,from,to,nStrips)
%function r=nma_trapezoidal(f,from,to,nStrips)
%
% integrates a function using trapezoidal rule using
% specific number of strips.
%
% INPUT:
% func : a string that repesents the function itself
% for example 'x*sin(x)'. The independent variable
% used in the string must be 'x' and no other letter.
%
% from: lower limit
% to : upper limit
% nStrips: number of strips to use
%
% OUTPUT
% I : The integral.
%
% Author: Nasser Abbasi
% May 3, 2003
I=0;
if(nStrips<=0)
error('number of strips must be > 0');
end
nPoints=nStrips+1;
X=linspace(from,to,nPoints);
h=X(2)-X(1);
for(i=1:length(X))
x=X(i);
f(i)=eval(func);
if(i==1 | i==length(X) )
I=I+f(i);
else
I=I+2*f(i);
end
end
I=I/2;
I=I*h;