Скачать Numerical Integration by the Romberg Method

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.