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.


function rect_left(integer rid, atom x, atom /*h*/)
    return rid(x)
end function
 
function rect_mid(integer rid, atom x, atom h)
    return rid(x+h/2)
end function
 
function rect_right(integer rid, atom x, atom h)
    return rid(x+h)
end function
 
function trapezium(integer rid, atom x, atom h)
    return (rid(x)+rid(x+h))/2
end function
 
function simpson(integer rid, atom x, atom h)
    return (rid(x)+4*rid(x+h/2)+rid(x+h))/6
end function
 
function cubed(atom x)
    return power(x,3)
end function
 
function recip(atom x)
    return 1/x
end function
 
function ident(atom x)
    return x
end function
 
function integrate(integer m_id, integer f_id, atom a, atom b, integer steps)
    atom accum = 0,
         h = (b-a)/steps 
    for i=0 to steps-1 do
        accum += m_id(f_id,a+h*i,h)
    end for
    return h*accum
end function
 
function smartp(atom N)
    if N=floor(N) then return sprintf("%d",N) end if
    string res = sprintf("%12f",round(N,1000000))
    if find('.',res) then
        res = trim_tail(res,"0")
        res = trim_tail(res,".")
    end if
    return res
end function
 
procedure test(sequence tests)
    string name
    atom a, b, steps, rid
    printf(1,"Function     Range     Iterations       L-Rect       M-Rect       R-Rect      Trapeze      Simpson\n")
    for i=1 to length(tests) do
        {name,a,b,steps,rid} = tests[i]
        printf(1,"  %-5s %6d - %-5d %10d  %12s %12s %12s %12s %12s\n",{name,a,b,steps,
                            smartp(integrate(rect_left, rid,a,b,steps)),
                            smartp(integrate(rect_mid,  rid,a,b,steps)),
                            smartp(integrate(rect_right,rid,a,b,steps)),
                            smartp(integrate(trapezium, rid,a,b,steps)),
                            smartp(integrate(simpson,   rid,a,b,steps))})
    end for
end procedure
 
constant tests = {{"x^3", 0,    1,     100, cubed},
                  {"1/x", 1,  100,    1000, recip},
                  {"x",   0, 5000, 5000000, ident},
                  {"x",   0, 6000, 6000000, ident}}
 
test(tests)
