Question 9.3: Romberg integration with comparison to the composite trapezo......
Romberg integration with comparison to the composite trapezoidal method.
Evaluate \int_{0}^{1} \frac{1}{(1+x)} d x using three levels of Romberg integration. Use an initial step size of h=1 (one subinterval). Compare your result with the exact answer. What number of subintervals would be required if you were to use the composite trapezoidal method to obtain the same level of accuracy?
Exact answer: The exact answer to this problem can be obtained analytically. The answer is: \ln (2)=0.69314718.
Script File
Romberg Integration
To carry out the numerical integration, a user-defined function, which is listed below, named Romberg is created.
Program 9-2: Function file. Romberg integration.
function IR = Romberg(Fun,a,b,Ni,Levels)
% Romberg numerically integrate using the Romberg integration method.
% Input Variables:
% Fun Name for the function to be integrated.
% (Fun is assumed to be written with element-by-element calculations.)
% a Lower limit of integration.
% b Upper limit of integration.
% Ni Initial number of subintervals.
% Levels Number of levels of Romberg integration.
% Output Variable:
% IR A matrix with the estimated values of the integral.
% Creating the first colwnn with the composite trapezoidal method:
for i = 1: Levels + 1
Nsubinter=Ni*2^ (i - 1);
IR(i,1)=trapezoidal(Fun,a,b,Nsubinter);
end
% Calculating the extrapolated values using Eq. (9.65):
for j = 2 : Levels + 1
for i=1: (Levels - j +2)
IR (i, j) = (4^ (j -1) *IR (i + 1, j -1) -IR (i, j -1)) / (4^ (j -1) -1) ;
end
end
The function Romberg is next used in the Command Window to determine the value of the integral \int_{0}^{1} \frac{1}{(1+x)} d x. The initial number of subintervals (for the first estimate with the composite trapezoidal method) is 1 .
>> format long
>> integ =@ (x) 1./(1+x);
>> IntVal = Romberg(integ,0,1,1,3)
IntVal =
0.75000000000000 0.69444444444444 0.69317460317460 0.69314747764483
0.70833333333333 0.69325396825397 0.69314790148123 0
0.69702380952381 0.69315453065453 0 0
0.69412185037185 0 0 0
The results show that with the composite trapezoidal method (first column) the most accurate value that is obtained (using eight subintervals) is accurate to two decimal places. The first-level Romberg integration (second column) increases the accuracy to four decimal places. The second-level Romberg integration (third column) increases the accuracy to six decimal places. The result from the third-level Romberg integration (fourth column) is also accurate to six decimal places, but the value is closer to the exact answer.
The number of calculations that was executed is 10 (4 in the composite trapezoidal method and 6 in the Romberg integration procedure). To obtain an estimate for the integral with an accuracy of six decimal places by applying only the composite trapezoidal method, the method has to be applied with 276 subintervals. This is shown below where trapezoidal is used in the Command Window.
>> ITrap = trapezoidal(integ,0,1,277)
ITrap =
0.69314799511374