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.


 
while 1
    read x$
    if x$ ="end" then print "**Over**": end
 
    read a, b, N, knownValue
 
    print " Function y ="; x$; " from "; a; " to "; b; " in "; N; " steps"
    print " Known exact value ="; knownValue
 
    areaLR = IntegralByLeftRectangle(   x$, a, b, N)
    areaRR = IntegralByRightRectangle(  x$, a, b, N)
    areaMR = IntegralByMiddleRectangle( x$, a, b, N)
    areaTr = IntegralByTrapezium(       x$, a, b, N)
    areaSi = IntegralBySimpsonRule(     x$, a, b, N)
 
    print "Left rectangle method   "; using( "##########.##########", areaLR); " diff "; knownValue-areaLR; tab(70); (knownValue-areaLR)/knownValue*100;" %"
    print "Right rectangle method  "; using( "##########.##########", areaRR); " diff "; knownValue-areaRR; tab(70); (knownValue-areaRR)/knownValue*100;" %"
    print "Middle rectangle method "; using( "##########.##########", areaMR); " diff "; knownValue-areaMR; tab(70); (knownValue-areaMR)/knownValue*100;" %"
    print "Trapezium  method       "; using( "##########.##########", areaTr); " diff "; knownValue-areaTr; tab(70); (knownValue-areaTr)/knownValue*100;" %"
    print "Simpson's Rule          "; using( "##########.##########", areaSi); " diff "; knownValue-areaSi; tab(70); (knownValue-areaSi)/knownValue*100;" %"
 
    print
 
wend
 
end
 
'------------------------------------------------------
    'we have N sizes, that gives us N+1 points
    'point 0 is a
    'point N is b
    'point i is xi =a +i *h
    'Often, precision is (sharper?) then single step area
    'So there should be EXACT number of steps, hence loop by integer i.
 
function IntegralByLeftRectangle( x$, a, b, N)
    h = ( b -a) /N
    s = 0
    for i = 0 to N -1
        x = a +i *h
        s = s + h *eval( x$)
    next
    IntegralByLeftRectangle = s
end function
 
function IntegralByRightRectangle( x$, a, b, N)
    h =( b -a) /N
    s = 0
    for i =1 to N
        x = a +i *h
        s = s + h *eval( x$)
    next
    IntegralByRightRectangle = s
end function
 
function IntegralByMiddleRectangle( x$, a, b, N)
    h =( b -a) /N
    s = 0
    for i =0 to N -1
        x = a +i *h +h /2
        s = s + h *eval( x$)
    next
    IntegralByMiddleRectangle = s
end function
 
function IntegralByTrapezium( x$, a, b, N)
'Formula is h*((f(a)+f(b))/2 + sum_{i=1}^{N-1} (f(x_i)))
    h  =( b -a) /N
    x  = a
    fa =eval( x$)
    x  =b
    fb =eval( x$)
    s = h *( fa +fb) /2
    for i =1 to N -1
        x = a +i *h
        s = s + h *eval( x$)
    next
    IntegralByTrapezium = s
end function
 
function IntegralBySimpsonRule( x$, a, b, N)
    'Simpson
    'N should be even.
    if N mod 2 then N =N +1
    'It really doesn't look right to double number of points from N to 2N -
    ' - this method is most accurate of all presented!
    'So we use NN as N/2, and N will be 2NN
    'Formula is h/6*( f(a)+f(b) + 4*(f(x_1)+f(x_3)+...+f(x_{2NN-1})+ 2*(f(x_2)+f(x_4)+...+f(x_{2NN-2})) )
    'Somehow I messed up h/6, h/3 and what is h, regarding "n=number of double intervals of size 2h"
    NN =N /2
 
    h  =( b -a) /N
    x  =a
    fa =eval (x$)
    x  =b
    fb =eval( x$)
    s = h /3 *( fa +fb)
    for i =1 to 2 *NN -1 step 2
        x = a +i *h
        s = s + h /3 *4 *eval( x$)  'odd points
    next
    for i =2 to 2 *NN -2 step 2
        x = a +i *h
        s = s + h /3 *2 *eval( x$)  'even points
    next
 
    IntegralBySimpsonRule = s
end function
 
'=======================================================
data "x^3",  0,    1,     100,          0.25
data "x^-1", 1,  100,    1000,          4.605170
data "x",    0, 5000,    1000,   12500000.0   '   should use 5 000 000 steps
data "x",    0, 6000,    1000,   18000000.0   '   should use 6 000 000 steps
data "end"
 
end
 