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 int_leftrect(sequence bounds, integer n, integer func_id)
    atom h, sum
    h = (bounds[2]-bounds[1])/n
    sum = 0
    for x = bounds[1] to bounds[2]-h by h do
        sum += call_func(func_id, {x})
    end for
    return h*sum
end function
 
function int_rightrect(sequence bounds, integer n, integer func_id)
    atom h, sum
    h = (bounds[2]-bounds[1])/n
    sum = 0
    for x = bounds[1] to bounds[2]-h by h do
        sum += call_func(func_id, {x+h})
    end for
    return h*sum
end function
 
function int_midrect(sequence bounds, integer n, integer func_id)
    atom h, sum
    h = (bounds[2]-bounds[1])/n
    sum = 0
    for x = bounds[1] to bounds[2]-h by h do
        sum += call_func(func_id, {x+h/2})
    end for
    return h*sum
end function
 
function int_trapezium(sequence bounds, integer n, integer func_id)
    atom h, sum
    h = (bounds[2]-bounds[1])/n
    sum = call_func(func_id, {bounds[1]}) + call_func(func_id, {bounds[2]})
    for x = bounds[1] to bounds[2]-h by h do
        sum += 2*call_func(func_id, {x})
    end for
    return h * sum / 2
end function
 
function int_simpson(sequence bounds, integer n, integer func_id)
    atom h, sum1, sum2
    h = (bounds[2]-bounds[1])/n
    sum1 = call_func(func_id, {bounds[1] + h/2})
    sum2 = 0
    for i = 1 to n-1 do
        sum1 += call_func(func_id, {bounds[1] + h * i + h / 2})
        sum2 += call_func(func_id, {bounds[1] + h * i})
    end for
    return h/6 * (call_func(func_id, {bounds[1]}) +
                  call_func(func_id, {bounds[2]}) + 4*sum1 + 2*sum2)
end function
 
function xp2d2(atom x)
    return x*x/2
end function
 
function logx(atom x)
    return log(x)
end function
 
function x(atom x)
    return x
end function
 
? int_leftrect({-1,1},1000,routine_id("xp2d2"))
? int_rightrect({-1,1},1000,routine_id("xp2d2"))
? int_midrect({-1,1},1000,routine_id("xp2d2"))
? int_simpson({-1,1},1000,routine_id("xp2d2"))
puts(1,'\n')
? int_leftrect({1,2},1000,routine_id("logx"))
? int_rightrect({1,2},1000,routine_id("logx"))
? int_midrect({1,2},1000,routine_id("logx"))
? int_simpson({1,2},1000,routine_id("logx"))
puts(1,'\n')
? int_leftrect({0,10},1000,routine_id("x"))
? int_rightrect({0,10},1000,routine_id("x"))
? int_midrect({0,10},1000,routine_id("x"))
? int_simpson({0,10},1000,routine_id("x"))