Simpsons Rule
Volume Number: 9
Issue Number: 10
Column Tag: Pascal workshop
Simpson’s Rule 
An ingenious method for approximating integrals
By Marek Hajek, Incline Village, Nevada
Note: Source code files accompanying article are located on MacTech CD-ROM orsource code disks.
About the author
Marek Hajek has been programming the Macintosh since 1989. He programmed
two and a half years for Sierra Software Innovations where he wrote several in-house
MacApp applications, participated in the development of SuperTEView, and the
relational database engine - Inside Out II. Currently, he is receiving his bachelor's
degree in Computer Science at the University of Nevada, Reno. He supports his college
education by making useful programming tools - sorting/ searching algorithms, and
custom development. He welcomes your comments on this article either by phone at
(702) 673-3341 or write to P.O. Box 7542, Incline Village, NV 89450.
Simpson's Rule, named after the great English mathematician Thomas Simpson, is
an ingenious method for approximating integrals. If you don't know what integrals are
used for, don't feel bad. Many college students who complete three semesters of
calculus may be able to “compute” an integral, but won't know its practical
application either. Computation of integrals is difficult to learn and easy to forget.
[Many years out of school, I can attest to this! - Ed.]
Integrals are essential to the modern world. Practical applications of the integral
are found in business, hydrostatics, highway construction, travel to the moon, solving
of differential equations, and other branches of science. The first computer ever built
was constructed to speed up ballistic missile trajectory calculations which meant
solving a lot of integrals. [Given the forces imposed on a missile (gravity, thrust, wind
resistance, etc.), integration is necessary to determine its path. - Tech. Ed.]
To illustrate the use of integrals, look at the curve in Figure 1.1a. The curve is
described by the equation (1+X4). I want to compute the area of the shaded region.
Notice, the area is between the curve, the x-axis, and the x-coordinates [-1,1]. The
integral that will compute the area of the shaded region is in figure 1.1b. Anybody
familiar with integrals will tell you that there is no known way to solve the integral
abstractly (the quick and easy way).
Figure 1.1a
Figure 1.1b
If an integral can be solved on an abstract level, the computation is relatively
easy. In practical applications, however, an integral can seldom be solved abstractly.
That's where the Simpson's Rule finds its use. Thomas Simpson invented an equation,
today called Simpson's Rule, which can be used to approximate an integral.
APPROXIMATION
The following line shows this equation in abstract form.
Looks complicated? First, take a look at figure 1.2.
Figure 1.2
In the approximation equation, the variables a and b are the boundaries of the
integral and correspond to -1 and 1 in figure 1.1a. The variable n is the number of
times the region under the curve is partitioned into smaller regions. You only have to
know two things about n. First, the larger n is, the more accurate the approximation.
And second, n must be a positive even integer (+2, 4, 6, ...). The function f(x) is the
function you are integrating. In my example, it is (1+X4). Whenever you encounter
f(x) in the equation, pass it the appropriate parameter. The parameters are the
x-coordinates of the partitions (X0, X1, X2, X3, , Xn).
[Simpson’s rule approximates the function on each subinterval of the partition
by a parabola that passes through the endpoints and the midpoint. The area under a
parabola is easily calculated. Adding up these areas gives an estimate of the integral. -
Tech. Ed.]
EXAMPLE COMPUTATION
To compute the integral in figure 1.1b, given four partitions (n=4), the
approximation looks like this:
Simplified:
Simplified:
Result: 2.1791. . .
SIMPSON'S RULE - PASCAL
Figure 1.3
I translated Simpson's Rule into several pascal functions. To help you see what
each function does, the equation is divided into three parts - Head, Twos/Fours, and
First/Last (Figure 1.3). Have fun!
CODE LISTING
{--------------------Main Program----------------------------}
PROGRAM Simpson;
(* Author - Marek Hajek *)
(* P.O. Box 7542 *)
(* Incline Village, NV 89450 *)
(* This program was written with Think Pascal 4.0.1 *)
USES
(* Make sure you include the sane library *)
Auxiliary, Sane;
CONST
kLowerLimit = -1; (* Corresponds to "a" *)
kUpperLimit = 1; (* Corresponds to "b" *)
kPartitions = 4; (* Corresponds to n = 4 *)
VAR
result: Extended;
(* The approximated result of the Integral *)
BEGIN
ShowText; (* Brings up the Think Pascal text window *)
result := ComputeIntegral(kLowerLimit, kUpperLimit,
kPartitions, IntegrandFunction);
writeln('Integral with lower/upper limits ', kLowerLimit : 0,
'/', kUpperLimit : 0, ', subintervals ', kPartitions : 0,
' is: ', result);
readln; (* Stop here before the text window disappears *)
END.
{--------------------ComputeIntegral-------------------------}
FUNCTION ComputeIntegral (lowerLimit, upperLimit: Extended;
partitionCount: LongInt;
FUNCTION IntegrandFunction (partitionCoordinate:
Extended): Extended): Extended;
(* The function ComputeIntegral calls the necessary *)
(* functions to compute the individual parts.*)
(* It returns the approximate result. *)
VAR
result: Extended;
head: Extended;
partitionIncrement: Extended;
partitionCoordinate: Extended;
index: LongInt;
BEGIN
head := ComputeHead(lowerLimit, upperLimit, partitionCount);
result := FirstAndLast(lowerLimit, upperLimit,
IntegrandFunction);
partitionIncrement :=
(upperLimit - lowerLimit) / partitionCount;
partitionCoordinate := lowerLimit;
(* The FOR loop computes the second part of the *)
(* integral -> Twos/Fours *)
FOR index := 1 TO partitionCount - 1 DO
BEGIN
(* Partition coordinate corresponds to X0, X1, X2,.....Xn *)
partitionCoordinate :=
partitionCoordinate + partitionIncrement;
(* Odd index means compute 4* f(x), even index *)
(* means compute 2 * f(x) *)
IF Odd(index) THEN
result := result +
4 * IntegrandFunction(partitionCoordinate)
ELSE
result := result +
2 * IntegrandFunction(partitionCoordinate)
END; (* FOR ... *)
ComputeIntegral := head * result;
END;
{------------------IntegrandFunction-------------------------}
FUNCTION IntegrandFunction (partitionCoordinate:
Extended): Extended;
(* The Integrand function is the function inside the *)
(* integral and needs to be defined by you. In my example, *)
(* the integrand function is (1+X4) and is translated *)
(* into pascal. The function takes one argument which is *)
(* the x coordinate of the partition. *)
BEGIN
{ This functions computes -> (X * X * X * X +1) }
IntegrandFunction :=
SQRT(XpwrI(partitionCoordinate, 4) + 1);
END;
{---------------------ComputeHead----------------------------}
FUNCTION ComputeHead (lowerLimit, upperLimit: Extended;
partitionCount: LongInt): Extended;
(* Computes the first part of the integral equation, *)
(* the Head. Corresponds to (b - a)/(3*n) *)
BEGIN
ComputeHead :=
(upperLimit - lowerLimit) / (3 * partitionCount);
END;
{----------------------FirstAndLast--------------------------}
FUNCTION FirstAndLast (lowerLimit, upperLimit: Extended;
FUNCTION IntegrandFunction (partitionCoordinate:
Extended): Extended): Extended;
(* Computes the third part of the integral, the *)
(* FIRST/LAST. Corresponds to [f(X0) + f(Xn) *)
BEGIN
FirstAndLast := IntegrandFunction(lowerLimit) +
IntegrandFunction(upperLimit);
END;