! ! =====> 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

#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) ); }

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*