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.


begin % compare some numeric integration methods                             %
 
    long real procedure leftRect ( long real procedure f
                                 ; long real value     a, b
                                 ; integer   value     n
                                 ) ;
    begin
        long real h, sum, x;
        h   := (b - a) / n;
        sum := 0;
        x   := a;
        while x <= b - h do begin
            sum := sum + (h * f(x));
            x   := x + h
        end;
        sum
    end leftRect ; 
 
    long real procedure rightRect ( long real procedure f
                                  ; long real value     a, b
                                  ; integer   value     n
                                  ) ;
    begin
        long real h, sum, x;
        h   := (b - a) / n;
        sum := 0;
        x   := a + h;
        while x <= b do begin
            sum := sum + (h * f(x));
            x   := x + h
        end;
        sum
    end rightRect ;
 
    long real procedure midRect ( long real procedure f
                                ; long real value     a, b
                                ; integer   value     n
                                ) ;
    begin
        long real h, sum, x;
        h   := (b - a) / n;
        sum := 0;
        x   := a;
        while x <= b - h do begin
            sum := sum + h * f(x + h / 2);
            x   := x + h
        end;
        sum
    end midRect ;
 
    long real procedure trapezium ( long real procedure f
                                  ; long real value     a, b
                                  ; integer   value     n
                                  ) ;
    begin
        long real h, sum, x;
        h   := (b - a) / n;
        sum := f(a) + f(b);
        x   := 1;
        while x <= n - 1 do begin
            sum := sum + 2 * f(a + x * h );
            x   := x + 1
        end;
        (b - a) / (2 * n) * sum
    end trapezium ;
 
    long real procedure simpson ( long real procedure f
                                ; long real value     a, b
                                ; integer   value     n
                                ) ;
    begin
        long real h, sum1, sum2, x;
        integer   limit;
        h     := (b - a) / n;
        sum1  := 0;
        sum2  := 0;
        limit := n - 1;
        for i := 0 until limit do sum1 := sum1 + f(a + h * i + h / 2);
        for i := 1 until limit do sum2 := sum2 + f(a + h * i);
        h / 6 * (f(a) + f(b) + 4 * sum1 + 2 * sum2)
    end simpson ;
 
    % tests the above procedures                                             %
    procedure testIntegrators1 ( string(3) value     legend
                               ; long real procedure f
                               ; long real value     lowerLimit
                               ; long real value     upperLimit
                               ; integer   value     iterations
                               ) ;
        write( r_format := "A", r_w := 20, r_d := 6, s_w := 0,
             , legend
             , leftRect(  f, lowerLimit, upperLimit, iterations )
             , rightRect( f, lowerLimit, upperLimit, iterations )
             , midRect(   f, lowerLimit, upperLimit, iterations )
             , trapezium( f, lowerLimit, upperLimit, iterations )
             , simpson(   f, lowerLimit, upperLimit, iterations )
             );
    procedure testIntegrators2 ( string(3) value     legend
                               ; long real procedure f
                               ; long real value     lowerLimit
                               ; long real value     upperLimit
                               ; integer   value     iterations
                               ) ;
        write( r_format := "A", r_w := 16, r_d := 2, s_w := 0,
             , legend
             , leftRect(  f, lowerLimit, upperLimit, iterations ), "    "
             , rightRect( f, lowerLimit, upperLimit, iterations ), "    "
             , midRect(   f, lowerLimit, upperLimit, iterations ), "    "
             , trapezium( f, lowerLimit, upperLimit, iterations ), "    "
             , simpson(   f, lowerLimit, upperLimit, iterations ), "    "
             );
 
    begin % task test cases                                                  %
        long real procedure xCubed   ( long real value x ) ; x * x * x;
        long real procedure oneOverX ( long real value x ) ; 1 / x;
        long real procedure xValue   ( long real value x ) ; x;
        write( "   "
             , "           left rect"
             , "          right rect"
             , "            mid rect"
             , "           trapezium"
             , "             simpson"
            );
        testIntegrators1( "x^3", xCubed,   0,    1,     100 );
        testIntegrators1( "1/x", oneOverX, 1,  100,    1000 );
        testIntegrators2( "x  ", xValue,   0, 5000, 5000000 );
        testIntegrators2( "x  ", xValue,   0, 6000, 6000000 )
    end
end.