Write functions to calculate the definite integral of a function  ƒ(x)  using all five of the following methods:

 rectangular
 left
 right
 midpoint
 trapezium
 Simpson's
 composite
Your functions should take in the upper and lower bounds (a and b), and the number of approximations to make in that range (n).

Assume that your example already has a function that gives values for  ƒ(x) .

Simpson's method is defined by the following pseudo-code:


 Pseudocode: Simpson's method, composite



procedure quad_simpson_composite(f, a, b, n)
    h := (b - a) / n
    sum1 := f(a + h/2)
    sum2 := 0

    loop on i from 1 to (n - 1)
        sum1 := sum1 + f(a + h * i + h/2)
        sum2 := sum2 + f(a + h * i)
                                            
    answer := (h / 6) * (f(a) + f(b) + 4*sum1 + 2*sum2)




Demonstrate your function by showing the results for:

   ƒ(x) = x3,       where   x   is     [0,1],       with           100 approximations.   The exact result is     0.25               (or 1/4)
   ƒ(x) = 1/x,     where   x   is   [1,100],     with        1,000 approximations.   The exact result is     4.605170+     (natural log of 100)
   ƒ(x) = x,         where   x   is   [0,5000],   with 5,000,000 approximations.   The exact result is   12,500,000
   ƒ(x) = x,         where   x   is   [0,6000],   with 6,000,000 approximations.   The exact result is   18,000,000


See also

   Active object for integrating a function of real time.
   Special:PrefixIndex/Numerical integration for other integration methods.


MODE F = PROC(LONG REAL)LONG REAL;
 
###############
## left rect ##
###############
 
PROC left rect = (F f, LONG REAL a, b, INT n) LONG REAL:
BEGIN
   LONG REAL h= (b - a) / n;
   LONG REAL sum:= 0;
   LONG REAL x:= a;
   WHILE x <= b - h DO
      sum := sum + (h * f(x));
      x +:= h
   OD;
   sum
END # left rect #; 
 
#################
## right rect  ##
#################
 
PROC right rect = (F f, LONG REAL a, b, INT n) LONG REAL:
BEGIN
   LONG REAL h= (b - a) / n;
   LONG REAL sum:= 0;
   LONG REAL x:= a + h;
   WHILE x <= b DO
      sum := sum + (h * f(x));
      x +:= h
   OD;
   sum
END # right rect #;
 
###############
## mid rect  ##
###############
 
PROC mid rect = (F f, LONG REAL a, b, INT n) LONG REAL:
BEGIN
   LONG REAL h= (b - a) / n;
   LONG REAL sum:= 0;
   LONG REAL x:= a;
   WHILE x <= b - h DO
      sum := sum + h * f(x + h / 2);
      x +:= h
   OD;
   sum
END # mid rect #;
 
###############
## trapezium ##
############### 
 
PROC trapezium = (F f, LONG REAL a, b, INT n) LONG REAL:
BEGIN
   LONG REAL h= (b - a) / n;
   LONG REAL sum:= f(a) + f(b);
   LONG REAL x:= 1;
   WHILE x <= n - 1 DO
      sum := sum + 2 * f(a + x * h );
      x +:= 1
   OD;
   (b - a) / (2 * n) * sum
END # trapezium #; 
 
#############
## simpson ##
############# 
 
PROC simpson = (F f, LONG REAL a, b, INT n) LONG REAL:
BEGIN
   LONG REAL h= (b - a) / n;
   LONG REAL sum1:= 0;
   LONG REAL sum2:= 0;
   INT limit:= n - 1;
   FOR i FROM 0 TO limit DO
      sum1 := sum1 + f(a + h * LONG REAL(i) + h / 2)
   OD;
   FOR i FROM 1 TO limit DO
      sum2 +:= f(a + h * LONG REAL(i))
   OD;
   h / 6 * (f(a) + f(b) + 4 * sum1 + 2 * sum2)
END # simpson #;
 
# test the above procedures #
PROC test integrators = ( STRING     legend
                        , F          function
                        , LONG REAL  lower limit
                        , LONG REAL  upper limit
                        , INT        iterations
                        ) VOID:
BEGIN
    print( ( legend
           , fixed(  left rect( function, lower limit, upper limit, iterations ), -20, 6 )
           , fixed( right rect( function, lower limit, upper limit, iterations ), -20, 6 )
           , fixed(   mid rect( function, lower limit, upper limit, iterations ), -20, 6 )
           , fixed(  trapezium( function, lower limit, upper limit, iterations ), -20, 6 )
           , fixed(    simpson( function, lower limit, upper limit, iterations ), -20, 6 )
           , newline
           )
         )
END; # test integrators #
print( ( "   "
       , "           left rect"
       , "          right rect"
       , "            mid rect"
       , "           trapezium"
       , "             simpson"
       , newline
       )
     );
test integrators( "x^3", ( LONG REAL x )LONG REAL: x * x * x, 0,     1,       100 );
test integrators( "1/x", ( LONG REAL x )LONG REAL: 1 / x,     1,   100,     1 000 );
test integrators( "x  ", ( LONG REAL x )LONG REAL: x,         0, 5 000, 5 000 000 );
test integrators( "x  ", ( LONG REAL x )LONG REAL: x,         0, 6 000, 6 000 000 );
 
SKIP