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: '
Print *, x0, x1
PRINT *

PRINT *, 'Number of panels (even): '
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: ');

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

REPEAT
Write('Enter number of panels for integration (even): ');
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]): ');