home *** CD-ROM | disk | FTP | other *** search
/ CD Actual 9 / CDACTUAL9.iso / share / Dos / VARIOS / pascal / MATH.SWG / 0002_Simpson's Approximation.pas < prev   
Encoding:
Pascal/Delphi Source File  |  1996-02-21  |  2.0 KB  |  84 lines
  1. { This simple Pascal program computes a definite integral using Simpson's
  2.   approximation.  Efficiency has been sacrificed for clarity, so it should
  3.   be easy to follow the logic of the program.                        - DJP }
  4.  
  5. Var                     { Declare variables   }
  6.   a : Integer;          {   Left point        }
  7.   b : Integer;          {   Right point       }
  8.   d : Real;             {   Delta x           }
  9.   i : Integer;          {   Iteration         }
  10.   n : Integer;          {   No. of intervals  }
  11.  
  12.   SubTotal : Real;
  13.  
  14. Function y(x:real):Real;   { Define your function here }
  15.   Begin
  16.   y:=1/x
  17.   End;
  18.  
  19. Function Coefficient(i:Integer):Integer;
  20.   Begin
  21.   If (i=0) or (i=n) Then
  22.     Coefficient:=1
  23.   Else
  24.     Coefficient:=(i Mod 2)*2+2
  25.  
  26.   { Notes:
  27.  
  28.     The MOD operater returns the remainder of a division.  This allows
  29.     us to determine if the partition is odd or even.
  30.  
  31.       <even> MOD 2 = 0
  32.       <odd>  MOD 2 = 1
  33.  
  34.     An examination of the coefficients of a typical approximation sum shows
  35.     an interesting pattern: Odd partitions have 4 as a coefficient and even
  36.     partitions have 2 as a coefficient.  The first and last partitions are
  37.     exceptions to this rule.  This pattern is used as a basis for
  38.     calculating the coefficient of a given partition.  }
  39.  
  40.   End;
  41.  
  42. Function xi(i:Integer):Real;
  43.   Begin
  44.   xi:=a+i*d
  45.   End;
  46.  
  47. Begin
  48.  
  49. a:=1;
  50. b:=2;
  51.  
  52. Repeat
  53.   Write('Subintervals? ');
  54.   ReadLn(n);
  55. Until (n Mod 2)=0; { Even number required }
  56.  
  57. d:=(b-1)/n;
  58.  
  59. WriteLn;
  60. WriteLn('  n      xi   f(xi)   c  cf(xi)');
  61. WriteLn('-------------------------------');
  62.  
  63. For i:=0 to n Do
  64.   Begin
  65.   WriteLn(i:3, xi(i):8:3, y(xi(i)):8:3, Coefficient(i):4,
  66.           Coefficient(i)*y(xi(i)):8:3);
  67.   SubTotal:=SubTotal+Coefficient(i)*y(xi(i))
  68.   End;
  69.  
  70. WriteLn;
  71. WriteLn('SubTotal',SubTotal:23:3);
  72. WriteLn('Result = ', (d/3)*SubTotal:0:50)
  73.  
  74. End.
  75.  
  76. { Quick Optimizations:
  77.  
  78.   a. Remove the WriteLn statement in the For loop.
  79.   b. Turn on 80x87 support.
  80.   c. Consolidate some of the procedures.
  81.  
  82.   peace/dp }
  83.  
  84.