1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
|
//Eg-11.9
//pg-488
clear
clc
a = 0;
b = 2;
h = b-a;
deff('out = func(in)','out = exp(-in^2)')
//From equations [30],[31],[32] & [33]
//Please note that the subscripts(i&j) we use here are different from that used in
//text book i.e they are increased by 1, because we cant give the index zero in //scilab. Therefore,
//Note : To get the results as in the book the lower limit of integration should be zero instead of 1 as in the book.
imax = 6;
jmax = 6;
I(1,1) = h/2*(func(a) + func(b));
I(2,1) = 1/2*(I(1,1) + h*func(a+h/2));
I(3,1) = 1/2*(I(2,1) + h/2*(func(a+h/4) + func(a+3*h/4)));
//From equation [33]
sum1 = 0;
for(j = 1:2:(2^3-1)) //Since we have to consider the odd terms only.
sum1 = sum1 + func(a+j*h/2^3);
end
I(4,1) = 1/2*(I(3,1) + h/2^2*sum1);
//Similarly
sum2 = 0;
for(j = 1:2:(2^4-1))
sum2 = sum2 + func(a+j*h/2^4);
end
I(5,1) = 1/2*(I(4,1) + h/2^3*sum2);
for(j = 2:5)
for(i = 1:imax-j)
I(i,j) = (4^(j-1)*I(i+1,j-1) - I(i,j-1))/(4^(j-1)-1);
end
end
printf(' The complete Romberg tableau is as follows\n')
disp(I)
printf('\n Therefore, the value of the integral is %f\n',I(1,5))
|
