[Back to MATH SWAG index] [Back to Main SWAG index] [Original]

    { This simple Pascal program computes a definite integral using Simpson's
      approximation.  Efficiency has been sacrificed for clarity, so it should
      be easy to follow the logic of the program.                        - DJP }

    Var                     { Declare variables   }
      a : Integer;          {   Left point        }
      b : Integer;          {   Right point       }
      d : Real;             {   Delta x           }
      i : Integer;          {   Iteration         }
      n : Integer;          {   No. of intervals  }

      SubTotal : Real;

    Function y(x:real):Real;   { Define your function here }
      Begin
      y:=1/x
      End;

    Function Coefficient(i:Integer):Integer;
      Begin
      If (i=0) or (i=n) Then
        Coefficient:=1
      Else
        Coefficient:=(i Mod 2)*2+2

      { Notes:

        The MOD operater returns the remainder of a division.  This allows
        us to determine if the partition is odd or even.

          <even> MOD 2 = 0
          <odd>  MOD 2 = 1

        An examination of the coefficients of a typical approximation sum shows
        an interesting pattern: Odd partitions have 4 as a coefficient and even
        partitions have 2 as a coefficient.  The first and last partitions are
        exceptions to this rule.  This pattern is used as a basis for
        calculating the coefficient of a given partition.  }

      End;

    Function xi(i:Integer):Real;
      Begin
      xi:=a+i*d
      End;

    Begin

    a:=1;
    b:=2;

    Repeat
      Write('Subintervals? ');
      ReadLn(n);
    Until (n Mod 2)=0; { Even number required }

    d:=(b-1)/n;

    WriteLn;
    WriteLn('  n      xi   f(xi)   c  cf(xi)');
    WriteLn('-------------------------------');

    For i:=0 to n Do
      Begin
      WriteLn(i:3, xi(i):8:3, y(xi(i)):8:3, Coefficient(i):4,
              Coefficient(i)*y(xi(i)):8:3);
      SubTotal:=SubTotal+Coefficient(i)*y(xi(i))
      End;

    WriteLn;
    WriteLn('SubTotal',SubTotal:23:3);
    WriteLn('Result = ', (d/3)*SubTotal:0:50)

    End.

    { Quick Optimizations:

      a. Remove the WriteLn statement in the For loop.
      b. Turn on 80x87 support.
      c. Consolidate some of the procedures.

      peace/dp }

[Back to MATH SWAG index] [Back to Main SWAG index] [Original]