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.


/*REXX pgm performs numerical integration using 5 different algorithms and show results.*/
numeric digits 20                                /*use twenty decimal digits precision. */
 
     do test=1  for 4;             say           /*perform the 4 different test suites. */
     if test==1  then do;    L= 0;     H=    1;     i=     100;     end
     if test==2  then do;    L= 1;     H=  100;     i=    1000;     end
     if test==3  then do;    L= 0;     H= 5000;     i= 5000000;     end
     if test==4  then do;    L= 0;     H= 6000;     i= 6000000;     end
     say center('test' test, 79, "═")            /*display a header for the test suite. */
     say '           left rectangular('L", "H', 'i")  ──► "         left_rect(L, H, i)
     say '       midpoint rectangular('L", "H', 'i")  ──► "     midpoint_rect(L, H, i)
     say '          right rectangular('L", "H', 'i")  ──► "        right_rect(L, H, i)
     say '                    Simpson('L", "H', 'i")  ──► "           Simpson(L, H, i)
     say '                  trapezium('L", "H', 'i")  ──► "         trapezium(L, H, i)
     end   /*test*/
exit                                             /*stick a fork in it,  we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
f:   parse arg y;  if test>2   then return y     /*choose the   "as─is"   function.     */
                   if test==1  then return y**3  /*   "    "     cube     function.     */
                                    return 1/y   /*   "    "  reciprocal     "          */
/*──────────────────────────────────────────────────────────────────────────────────────*/
left_rect:     procedure expose test; parse arg a,b,#;     $= 0;                h= (b-a)/#
                             do x=a      by h  for #;      $= $ + f(x)
                             end   /*x*/
               return $*h/1
/*──────────────────────────────────────────────────────────────────────────────────────*/
midpoint_rect: procedure expose test; parse arg a,b,#;     $= 0;                h= (b-a)/#
                             do x=a+h/2  by h  for #;      $= $ + f(x)
                             end   /*x*/
               return $*h/1
/*──────────────────────────────────────────────────────────────────────────────────────*/
right_rect:    procedure expose test; parse arg a,b,#;     $= 0;                h= (b-a)/#
                             do x=a+h    by h  for #;      $= $ + f(x)
                             end   /*x*/
               return $*h/1
/*──────────────────────────────────────────────────────────────────────────────────────*/
Simpson:       procedure expose test; parse arg a,b,#;                          h= (b-a)/#
               hh= h/2;                                    $= f(a + hh)
               @= 0;         do x=1  for #-1; hx=h*x + a;  @= @ + f(hx)
                                                           $= $ + f(hx + hh)
                             end   /*x*/
 
               return h * (f(a) + f(b) + 4*$ + 2*@)  /  6
/*──────────────────────────────────────────────────────────────────────────────────────*/
trapezium:     procedure expose test; parse arg a,b,#;     $= 0;                h= (b-a)/#
                             do x=a  by h  for #;          $= $ + (f(x) + f(x+h))
                             end   /*x*/
               return $*h/2