Computers in Engineering WWW Site - Example 17.1

Example 17.1


FORTRAN version

!
! =====> Program - P124.F90
!
        PROGRAM Intgrt
!  Program to demonstrate numerical integration.


!       Define function return types!!
        REAL            MidPt, Trapzd, Simpsn

        INTEGER         Panels
        REAL            x0, x1

        PRINT *, 'This is Program >> P124 = Integration Demo'
!
!     Tell program where data for  READ *  is coming from
      OPEN(UNIT=5, FILE='P124.DAT')      ! UNIT=5 is the default input
!
        PRINT *, 'Demonstration of Numerical Integration'
        PRINT *,'Function = 1/(1 + x*x)  for example'
        PRINT *
        PRINT *, 'Low limit, high limit: '
        READ *, x0, x1
        Print *, x0, x1
        PRINT *

        PRINT *, 'Number of panels (even): '
        READ *, Panels
        Print *, Panels
        PRINT *

        PRINT *, 'Midpoint:  ', MidPt(x0, x1, Panels)
        PRINT *, 'Trapezoid: ', Trapzd(x0, x1, Panels)
        PRINT *, 'Simpson''s  ', Simpsn(x0, x1, Panels)

        STOP
        END


        REAL FUNCTION F(x)
!       The function to be integrated.

        F = 1/(1 + x*x)
        RETURN
        END


        REAL FUNCTION MidPt(x0, x1, N)
!  Estimate the area under F between x0 and x1
!  using the Midpoint Rule with N panels.

        REAL            x0, x1,    &

                        x, h, Sum
        INTEGER         N


        h = (x1 - x0) / N
        Sum = 0.0

!       Add area of each panel, assuming unit width
!       and multiplying by the true width h later.
        x = x0+h/2.0   ! Initial start is middle of first panel
L1:     DO k= 1,N
          Sum = Sum + F(x)
          x = x + h
        END DO L1

        MidPt = h*Sum

        RETURN
        END


        REAL FUNCTION Trapzd(x0, x1, N)
!  Estimate the area under F between x0 and x1
!  using the Trapezoidal Rule with N panels.

        REAL            x0, x1,    &

                        x, h, Sum
        INTEGER         N


        h = (x1 - x0) / N
        Sum = (F(x0) + F(x1))/2.0

!       Add contributions from interior trapezoids.
        x = x0+h
L1:     DO k = 1,N-1
          Sum = Sum + F(x)
          x = x + h
        END DO L1

        Trapzd = h*Sum

        RETURN
        END


        REAL FUNCTION Simpsn(x0, x1, N)
!  Estimate the area under F between x0 and x1
!  using Simpson's Rule with N/2 panels.

        REAL            x0, x1,    &

                        x, h, Sum
        INTEGER         N


        h = (x1 - x0) / N
        Sum = F(x0) + F(x1)

!       Add 4 * odd-indexed f values and 2 * even-
!       indexed f values to Sum
        x = x0
L1:     DO I = 1, N-1
          x = x + h

          IF (Mod(I, 2) == 0) THEN
            Sum = Sum + 2.0*F(x)
          ELSE
            Sum = Sum + 4.0*F(x)
          END IF
        END DO L1

        Simpsn = h/3.0*Sum

        RETURN
        END

C version

#include <stdio.h>
#include <math.h>

/* Demonstration of numerical integration algorithms */

/*Program to demonstrate numerical integration methods.  A definite
  integral is an area under a curve, approximated by the sum of areas
  of small panels.


  The Midpoint, Trapezoidal, and Simpson's rules each use different
  panel shapes.  Each method has a similar computational cost, but
  Simpson's method has the edge in accuracy, being exact for
  quadratic functions (polynomials of degree 2).
*/



double f(double x){
/* The function to be integrated.  f(x)= x^3 + 5x -10  */

  return( pow(x,3) + 5.0*x - 10.0 );
}


double area_f(double x0, double x1){
/* Returns the definite integral of F between x0 and x1. */
  double temp1, temp2;
  temp1 = pow(x0,4)/4.0 + pow(x0,2)*5.0/2.0 - 10.0*x0;
  temp2 = pow(x1,4)/4.0 + pow(x1,2)*5.0/2.0 - 10.0*x1;
  return(temp2 - temp1);
}

void swap(double *a, double *b){
/* this function exchanges the values of two double variables */

  double temp;
  temp = *a;
  *a = *b;
  *b = temp;
}

double midpoint(double x0, double x1, int n){
/* Estimate the area under F between x0 and x1 using the Midpoint
  Rule with N panels.  The area is approximated by N rectangles of
  height determined by the value of F at the center of each panel. */

  double h, sum;


/* Calculate panel width h; initialize Sum. */
  h= (x1 - x0) / n;
  x0= x0 + h/2;
  sum= 0;

/* Add area of each panel, assuming unit width and multiplying
    by the true width h only once later. */
  do {
      sum= sum + f(x0);
      x0= x0 + h;
      }while(x0 < x1);


  return(sum * h);
}


double trapezoid(double x0, double x1, int n){
/* Estimate the area under F between x0 and x1 using the Trapezoidal
  Rule with N trapezoidal panels.  The height of the two parallel sides
  of each trapezoid is given by the value of F at either edge of
  the panel. */

  double h, sum;
  int i;


/* Calculate panel width and initialize sum. */
  h= (x1 - x0) / n;
  sum= f(x0)/2;


/* Remember that the two endpoints only contribute half as
   much to the final sum as the inner points do.  This is
   handled here by initializing the sum to f(x0)/2, and after
   adding up the contributions of all the other points, we must
   subtract half the contribution of the endpoint. */

  for(i=1; i<=n; i++)
      sum= sum + f(x0 + i*h);
  sum= sum - f(x1)/2;
  return(sum * h);
}



double simpson(double x0, double x1, int n){
/* Estimate the area under F between x0 and x1 using Simpson's Rule
  with N panels.  The area is approximated by that under N/2
  parabolic segments covering each of the N/2 pairs of adjacent
  panels. */

  double h, x, sum;
  int i;



/* Ensure even number of panels. If remainder after dividing by two is
   not zero, the number was odd, and should be incremented by 1 */

  if( (n%2) != 0 ) n++;

/* Calculate panel width h; initialize Sum. */
  h= (x1 - x0) / n;
  sum= f(x0);


/*  The function values at the endpoints contribute 1 each,
    the function values at the even points contribute 4 each,
    and the function values at the interior odd points; 2 each. */

  for(i=2; i<=n; i+=2){

    x= x0 + i*h;
    sum= sum + 4*f(x-h) + 2*f(x);
  }
  sum = sum - f(x1);
  return(h*sum / 3);
}



main(){

  double x0, x1, area, truearea;
  int i, panels;
  printf("Numerical Integration Demonstration\n\n");
  printf("  This program will approximate the definite integral of\n");
  printf("x^3 + 5x - 10 using three methods - the midpoint rule, the\n");
  printf("trapezoid rule and Simpson's rule.  You will be prompted for\n");
  printf("the beginning and end points as well as the number of panels.\n");
  printf("Note that the methods only vary in the shapes of the areas\n");
  printf("they use to make the approximations.\n\n");
  printf("Enter lower limit (x0): ");
  scanf("%lf",&x0);
  printf("Enter higher limit (x1): ");
  scanf("%lf",&x1);

  if(x0 > x1) swap(&x0, &x1);

  do{
     printf("Enter number of panels for integration (even): ");
     scanf("%d",&panels);
     } while( (panels % 2) != 0 && panels > 0);

  printf("\n\n Method    Panels    Area          True          Error\n");
  for(i=1; i<= 53; i++)
    printf("-");
  printf("\n");

  truearea= area_f(x0, x1);

    area= midpoint(x0, x1, panels);
    printf(" Midpoint  %6d    %12.4lf  %12.4lf  %12.4lf\n",panels,
           area, truearea, (area - truearea) );

    area= trapezoid(x0, x1, panels);
    printf(" Trapezoid %6d    %8.4lf  %8.4lf  %8.4lf\n",panels,
           area, truearea, (area - truearea) );

    area= simpson(x0, x1, panels);
    printf(" Simpson   %6d    %8.4lf  %8.4lf  %8.4lf\n",panels,
           area, truearea, (area - truearea) );
}

Pascal version

PROGRAM IntegrateDemo;
{ Program to demonstrate numerical integration methods.  A definite
  integral is an area under a curve, approximated by the sum of areas
  of small panels.


  The Midpoint, Trapezoidal, and Simpson's rules each use different
  panel shapes.  Each method has a similar computational cost, but
  Simpson's method has the edge in accuracy, being exact for
  quadratic functions (polynomials of degree 2).

  This program was compiled with Turbo Pascal V5.5, but should be
  compilable with Turbo Pascal versions 3.0 and above.  Deviations from
  standard Pascal are also minimal, particularly in the coding of the
  algorithms themselves. }



VAR
  Panels:       Integer;
  X0, X1,
  Area,
  TrueArea:     Real;
  Ch:           Char;



FUNCTION F(x: Real): Real;
{ The function to be integrated. }

BEGIN
  F:= 1/(1 + Sqr(x));
END;  { FUNCTION F }



FUNCTION Area_F(x0, x1: Real): Real;
{ Returns the definite integral of F between x0 and x1. }

BEGIN
  Area_F:= ArcTan(x1) - ArcTan(x0);
END;  { FUNCTION Area_F }



PROCEDURE Swap(VAR x1, x2: Real);
{ Exchanges two real values x1 and x2. }

VAR
  Temp: Real;

BEGIN
  Temp:= x1;
  x1:= x2;
  x2:= Temp;
END;  { PROCEDURE Swap }



FUNCTION MidPoint(x0, x1: Real; N: Integer): Real;
{ Estimate the area under F between x0 and x1 using the Midpoint
  Rule with N panels.  The area is approximated by N rectangles of
  height determined by the value of F at the center of each panel. }

VAR
  h,
  Sum:          Real;

BEGIN
  IF x0 > x1 THEN
    Swap(x0, x1);

  { Calculate panel width h; initialize Sum. }
  h:= (x1 - x0) / N;
  x0:= x0 + h/2;
  Sum:= 0;

  { Add area of each panel, assuming unit width and multiplying
    by the true width h only once later. }
  WHILE x0 < x1 DO
    BEGIN
      Sum:= Sum + F(x0);
      x0:= x0 + h;
    END;  { WHILE }

  MidPoint:= h*Sum;

END;  { FUNCTION MidPoint }



FUNCTION Trapezoid(x0, x1: Real; N: Integer): Real;
{ Estimate the area under F between x0 and x1 using the Trapezoidal
  Rule with N trapezoidal panels.  The height of the two parallel sides
  of each trapezoid is given by the value of F at either edge of
  the panel. }

VAR
  I:            Integer;
  h,
  Sum:          Real;

BEGIN
  IF x0 > x1 THEN
    Swap(x0, x1);

  { Calculate panel width and initialize sum. }
  h:= (x1 - x0) / N;
  Sum:= (F(x0) + F(x1))/2;

  { Add contributions from interior trapezoids. }
  FOR I:= 1 TO Pred(N) DO
    BEGIN
      x0:= x0 + h;
      Sum:= Sum + F(x0);
    END;  { FOR I }

  Trapezoid:= h*Sum;
END;  { FUNCTION Trapezoid }



FUNCTION Simpson(x0, x1: Real; N: Integer): Real;
{ Estimate the area under F between x0 and x1 using Simpson's Rule
  with N panels.  The area is approximated by that under N/2
  parabolic segments covering each of the N/2 pairs of adjacent
  panels. }

VAR
  I:            Integer;
  h,
  Sum:          Real;

BEGIN
  IF x0 > x1 THEN
    Swap(x0, x1);

  { Ensure even number of panels. }
  IF Odd(N) THEN
    N:= Succ(N);

  { Calculate panel width h; initialize Sum. }
  h:= (x1 - x0) / N;
  Sum:= F(x0) + F(x1);

  { If the points between x0 and x1 are indexed 0, 1, ..., N,
    then Simpson's Rule requires adding F0, FN (the interval extrema),
    four times points with odd subscripts (the midpoints of each 3-point
    double panel), and two times points with even subscripts (the
    endpoints of the interior double panels). }
  FOR I:= 1 TO Pred(N) DO
    BEGIN
      x0:= x0 + h;
      IF Odd(I) THEN
        Sum:= Sum + 4*F(x0)
      ELSE
        Sum:= Sum + 2*F(x0);
    END;  { FOR I }

  { Finish up by multiplying Sum by one third the panel width, h. }
  Simpson:= h/3*Sum;
END;  { FUNCTION Simpson }



BEGIN  { Main }
  WriteLn('Numerical Integration Demonstration');

  REPEAT
    WriteLn;
    Write('Enter low limit, high limit: ');
    ReadLn(x0, x1);

    IF x0 > x1 THEN
      Swap(x0, x1);

    REPEAT
      Write('Enter number of panels for integration (even): ');
      ReadLn(Panels);
    UNTIL (Panels > 0) AND NOT Odd(Panels);

    WriteLn; WriteLn;
    WriteLn('Method':15, 'Panels':8, 'Area':12, 'True':12, 'Error':12);

    TrueArea:= Area_F(x0, x1);

    WriteLn;
    Area:= MidPoint(x0, x1, Panels);
    WriteLn('Midpoint':15, Panels:8, Area:12:8, TrueArea:12:8,
        (Area - TrueArea):12:8);

    Area:= Trapezoid(x0, x1, Panels);
    WriteLn('Trapezoidal':15, Panels:8, Area:12:8, TrueArea:12:8,
        (Area - TrueArea):12:8);

    Area:= Simpson(x0, x1, Panels);
    WriteLn('Simpson':15, Panels:8, Area:12:8, TrueArea:12:8,
        (Area - TrueArea):12:8);

    WriteLn; WriteLn;
    Write('Another run? (y/[N]): ');
    ReadLn(Ch);
  UNTIL NOT ((Ch = 'Y') OR (Ch = 'y'));
END.

Last modified: 08/07/97