home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / TURBOPAS / PAS_ENG.ZIP / ROMB3.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1985-07-19  |  1.9 KB  |  95 lines
  1. program romb3;  { -> 287 }
  2. { integration by the romberg method }
  3.  
  4. const tol  = 1.0E-4;
  5. var done  : boolean;
  6.  sumt  : real;
  7.  sum,upper,lower : real;
  8.  
  9.  
  10. function fx(x: real): real;
  11. { find f(x)= 1/sqrt(x); watch out for x=0 }
  12. begin
  13.   fx:=1.0/sqrt(x)
  14. end;
  15.  
  16. procedure romb(
  17. lower,upper,tol: real;
  18.    var ans: real);
  19. { numerical integration by the Romberg method }
  20. var
  21.  nx   : array[1..16] of integer;
  22.  t   : array[1..136] of real;
  23.  done,error  : boolean;
  24.  pieces,nt,i,ii,n,nn,
  25.  l,ntra,k,m,j  : integer ;
  26.  delta_x,c,sum,fotom,x : real ;
  27. begin
  28.   done:=false;
  29.   error:=false;
  30.   pieces:=1;
  31.   nx[1]:=1;
  32.   delta_x:=(upper-lower)/pieces;
  33.   c:=(fx(lower)+fx(upper))*0.5;
  34.   t[1]:=delta_x*c;
  35.   n:=1;
  36.   nn:=2;
  37.   sum:=c;
  38.   repeat
  39.     n:=n+1;
  40.     fotom:=4.0;
  41.     nx[n]:=nn;
  42.     pieces:=pieces*2;
  43.     l:=pieces-1;
  44.     delta_x:=(upper-lower)/pieces;
  45.     { compute trapezoidal sum for 2^(n-1)+1 points }
  46.     for ii:=1 to (l+1) div 2 do
  47.       begin
  48.  i:=ii*2-1;
  49.  x:=lower+i*delta_x;
  50.  sum:=sum+fx(x)
  51.       end;
  52.     t[nn]:=delta_x*sum;
  53.     ntra:=nx[n-1];
  54.     k:=n-1;
  55.     { compute n-th row of T array }
  56.     for m:=1 to k do
  57.       begin
  58.  j:=nn+m;
  59.  nt:=nx[n-1]+m-1;
  60.  t[j]:=(fotom*t[j-1]-t[nt])/(fotom-1.0);
  61.  fotom:=fotom*4.0
  62.       end;
  63.     if n>4 then
  64.       begin
  65.  if t[nn+1]<>0.0 then
  66.    if (abs(t[ntra+1]-t[nn+1])<=abs(t[nn+1]*tol))
  67.      or (abs(t[nn-1]-t[j])<=abs(t[j]*tol)) then
  68.        done:=true
  69.    else if n>15 then
  70.      begin
  71.        done:=true;
  72.        error:=true
  73.      end
  74.  end; { if n>4 }
  75.     nn:=j+1
  76.   until done;
  77.   ans:=t[j]
  78. end;  { ROMBERG }
  79.  
  80. begin  { main program }
  81.   clrscr;
  82.   lower:=0.1;
  83.   upper:=1.0;
  84.   writeln;
  85.   sumt:=0.0;
  86.   writeln('new area total area lower upper limits');
  87.   repeat
  88.     romb(lower,upper,tol,sum);
  89.     upper:=lower;
  90.     lower:=0.1*upper;
  91.     sumt:=sumt+sum;
  92.     writeln(sum:9:6,'    ',sumt:9:5,'    ',lower,'    ',upper)
  93.   until abs(sum)<tol
  94. end.
  95.