Category : Pascal Source Code
Archive   : ENGR.ZIP
Filename : ROMB1.PAS

program romb1; { -> 281 }
{ integration by the romberg method }

const tol = 1.0E-4;
var done : boolean;
sum,upper,lower : real;

external procedure cls;

function fx(x: real): real;
{ find f(x)= 1.0/x; watch out for x=0 }
begin
fx:=1.0/x
end;

procedure romb(function f(x:real): real;
lower,upper,tol: real;
var ans: real);
{ numerical integration by the Romberg method }
var
nx : array[1..16] of integer;
t : array[1..136] of real;
done,error : boolean;
pieces,nt,i,ii,n,nn,
l,ntra,k,m,j : integer ;
delta_x,c,sum,fotom,x : real ;
begin
done:=false;
error:=false;
pieces:=1;
nx[1]:=1;
delta_x:=(upper-lower)/pieces;
c:=(f(lower)+f(upper))*0.5;
t[1]:=delta_x*c;
n:=1;
nn:=2;
sum:=c;
repeat
n:=n+1;
fotom:=4.0;
nx[n]:=nn;
pieces:=pieces*2;
l:=pieces-1;
delta_x:=(upper-lower)/pieces;
{ compute trapezoidal sum for 2^(n-1)+1 points }
for ii:=1 to (l+1) div 2 do
begin
i:=ii*2-1;
x:=lower+i*delta_x;
sum:=sum+f(x)
end;
t[nn]:=delta_x*sum;
write(pieces:5,t[nn]);
ntra:=nx[n-1];
k:=n-1;
{ compute n-th row of T array }
for m:=1 to k do
begin
j:=nn+m;
nt:=nx[n-1]+m-1;
t[j]:=(fotom*t[j-1]-t[nt])/(fotom-1.0);
fotom:=fotom*4.0
end;
writeln(j:4,t[j]);
if n>4 then
begin
if t[nn+1]<>0.0 then
if (abs(t[ntra+1]-t[nn+1])<=abs(t[nn+1]*tol))
or (abs(t[nn-1]-t[j])<=abs(t[j]*tol)) then
done:=true
else if n>15 then
begin
done:=true;
error:=true
end
end; { if n>4 }
nn:=j+1
until done;
ans:=t[j]
end; { ROMBERG }

begin { main program }
cls;
lower:=1.0;
upper:=9.0;
writeln;
romb(fx,lower,upper,tol,sum);
writeln;
writeln(chr(7),'Area= ',sum)
end.