home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / TURBOPAS / PAS_ENG.ZIP / ERFSIMP.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1985-07-19  |  1.6 KB  |  78 lines

  1. program erfsimp;  { -> 323 }
  2.  
  3. { integration by Simpson's method }
  4.  
  5. const tol  = 1.0E-4;
  6.  pi  = 3.141593;
  7.  
  8. var done  : boolean;
  9.  sum,upper,lower,
  10.  erf,twopi : real;
  11.  
  12.  
  13. function fx(x: real): real;
  14. begin
  15.   fx:=exp(-x*x)
  16. end;  { function fx }
  17.  
  18.  
  19. procedure simps(
  20.   lower,upper,tol : real;
  21.   var sum  : real);
  22.  
  23. { numerical integration by Simpson's rule }
  24. { function is f, limits are lower and upper }
  25. { with number of regions equal to pieces }
  26. { partition is delta_x, answer is sum }
  27.  
  28. var i  : integer;
  29.  x,delta_x,even_sum,
  30.  odd_sum,end_sum,
  31.  sum1  : real;
  32.  pieces  : integer;
  33. begin
  34.   pieces:=2;
  35.   delta_x:=(upper-lower)/pieces;
  36.   odd_sum:=fx(lower+delta_x);
  37.   even_sum:=0.0;
  38.   end_sum:=fx(lower)+fx(upper);
  39.   sum:=(end_sum+4.0*odd_sum)*delta_x/3.0;
  40.   writeln(pieces:5,sum);
  41.   repeat
  42.     pieces:=pieces*2;
  43.     sum1:=sum;
  44.     delta_x:=(upper-lower)/pieces;
  45.     even_sum:=even_sum+odd_sum;
  46.     odd_sum:=0.0;
  47.     for i:=1 to pieces div 2 do
  48.       begin
  49.  x:=lower+delta_x*(2.0*i-1.0);
  50.  odd_sum:=odd_sum+fx(x)
  51.       end;
  52.     sum:=(end_sum+4.0*odd_sum+2.0*even_sum)*delta_x/3.0;
  53.   until abs(sum-sum1)<=abs(tol*sum1)
  54. end; { simps }
  55.  
  56. begin  { main program }
  57.   clrscr;
  58.   done:=false;
  59.   twopi:=2.0/sqrt(pi);
  60.   lower:=0.0;
  61.  
  62.   repeat
  63.     writeln;
  64.     writeln('Erf? ');
  65.     readln(upper);
  66.     if upper<0.0 then done:=true
  67.     else if upper=0.0 then
  68.       writeln('Erf of 0.0 is 0.0')
  69.  
  70.  else { upper>0 }
  71.    begin
  72.      simps(lower,upper,tol,sum);
  73.      erf:=twopi*sum;
  74.      writeln('Erf of ',upper:7:2,', is ',erf:8:4)
  75.         end
  76.     until done
  77. end.
  78.