home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / TURBOPAS / PAS_ENG.ZIP / SIMP1.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1985-07-19  |  1.3 KB  |  63 lines
  1. program simp1;  { -> 273 }
  2. { integration by Simpson's method }
  3.  
  4. const tol  = 1.0E-4;
  5. var sum,upper,lower : real;
  6.  
  7.  
  8. function fx(x: real): real;
  9. { find f(x)=1/x }
  10. { watch out for x=0 }
  11.  
  12. begin
  13.   fx:=1.0/x
  14. end; { function fx }
  15.  
  16. procedure simps(
  17.   lower,upper,tol : real;
  18.   var sum  : real);
  19.  
  20. { numerical integration by Simpson's rule }
  21. { function is f (as paramater), limits are lower and upper }
  22. { with number of regions equal to pieces }
  23. { partition is delta_x, answer is sum }
  24.  
  25. var i  : integer;
  26.  x,delta_x,even_sum,
  27.  odd_sum,end_sum,
  28.  sum1  : real;
  29.  pieces  : integer;
  30. begin
  31.   pieces:=2;
  32.   delta_x:=(upper-lower)/pieces;
  33.   odd_sum:=fx(lower+delta_x);
  34.   even_sum:=0.0;
  35.   end_sum:=fx(lower)+fx(upper);
  36.   sum:=(end_sum+4.0*odd_sum)*delta_x/3.0;
  37.   writeln(pieces:5,sum);
  38.   repeat
  39.     pieces:=pieces*2;
  40.     sum1:=sum;
  41.     delta_x:=(upper-lower)/pieces;
  42.     even_sum:=even_sum+odd_sum;
  43.     odd_sum:=0.0;
  44.     for i:=1 to pieces div 2 do
  45.       begin
  46.  x:=lower+delta_x*(2.0*i-1.0);
  47.  odd_sum:=odd_sum+fx(x)
  48.       end;
  49.     sum:=(end_sum+4.0*odd_sum+2.0*even_sum)*delta_x/3.0;
  50.     writeln(pieces:5,sum)
  51.   until abs(sum-sum1)<=abs(tol*sum)
  52. end; { simps }
  53.  
  54. begin  { main program }
  55.   clrscr;
  56.   lower:=1.0;
  57.   upper:=9.0;
  58.   writeln;
  59.   simps(lower,upper,tol,sum);
  60.   writeln;
  61.   writeln(chr(7),'area= ',sum)
  62. end.
  63.