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.