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.


 
     1000 PRINT "F(X)";" FROM";"   TO";"       L-Rect";"       M-Rect";"       R-Rect ";"      Trapez";"      Simpson"
     1010 fromval:=0
     1020 toval:=1
     1030 PRINT "X^3 ";
     1040 PRINT USING "#####": fromval;
     1050 PRINT USING "#####": toval;
     1060 PRINT USING "###.#########": numint(f1, "L", fromval, toval, 100);
     1070 PRINT USING "###.#########": numint(f1, "R", fromval, toval, 100);
     1080 PRINT USING "###.#########": numint(f1, "M", fromval, toval, 100);
     1090 PRINT USING "###.#########": numint(f1, "T", fromval, toval, 100);
     1100 PRINT USING "###.#########": numint(f1, "S", fromval, toval, 100)
     1110 //
     1120 fromval:=1
     1130 toval:=100
     1140 PRINT "1/X ";
     1150 PRINT USING "#####": fromval;
     1160 PRINT USING "#####": toval;
     1170 PRINT USING "###.#########": numint(f2, "L", fromval, toval, 1000);
     1180 PRINT USING "###.#########": numint(f2, "R", fromval, toval, 1000);
     1190 PRINT USING "###.#########": numint(f2, "M", fromval, toval, 1000);
     1200 PRINT USING "###.#########": numint(f2, "T", fromval, toval, 1000);
     1210 PRINT USING "###.#########": numint(f2, "S", fromval, toval, 1000)
     1220 fromval:=0
     1230 toval:=5000
     1240 PRINT "X   ";
     1250 PRINT USING "#####": fromval;
     1260 PRINT USING "#####": toval;
     1270 PRINT USING "#########.###": numint(f3, "L", fromval, toval, 5000000);
     1280 PRINT USING "#########.###": numint(f3, "R", fromval, toval, 5000000);
     1290 PRINT USING "#########.###": numint(f3, "M", fromval, toval, 5000000);
     1300 PRINT USING "#########.###": numint(f3, "T", fromval, toval, 5000000);
     1310 PRINT USING "#########.###": numint(f3, "S", fromval, toval, 5000000)
     1320 //
     1330 fromval:=0
     1340 toval:=6000
     1350 PRINT "X   ";
     1360 PRINT USING "#####": fromval;
     1370 PRINT USING "#####": toval;
     1380 PRINT USING "#########.###": numint(f3, "L", fromval, toval, 6000000);
     1390 PRINT USING "#########.###": numint(f3, "R", fromval, toval, 6000000);
     1400 PRINT USING "#########.###": numint(f3, "M", fromval, toval, 6000000);
     1410 PRINT USING "#########.###": numint(f3, "T", fromval, toval, 6000000);
     1420 PRINT USING "#########.###": numint(f3, "S", fromval, toval, 6000000)
     1430 END
     1440 //
     1450 FUNC numint(FUNC f, type$, lbound, rbound, iters) CLOSED
     1460   delta:=(rbound-lbound)/iters
     1470   integral:=0
     1480   CASE type$ OF
     1490   WHEN "L", "T", "S"
     1500     actval:=lbound
     1510   WHEN "M"
     1520     actval:=lbound+delta/2
     1530   WHEN "R"
     1540     actval:=lbound+delta
     1550   OTHERWISE
     1560     actval:=lbound
     1570   ENDCASE
     1580   FOR n:=0 TO iters-1 DO
     1590     CASE type$ OF
     1600     WHEN "L", "M", "R"
     1610       integral:+f(actval+n*delta)*delta
     1620     WHEN "T"
     1630       integral:+delta*(f(actval+n*delta)+f(actval+(n+1)*delta))/2
     1640     WHEN "S"
     1650       IF n=0 THEN
     1660         sum1:=f(lbound+delta/2)
     1670         sum2:=0
     1680       ELSE
     1690         sum1:+f(actval+n*delta+delta/2)
     1700         sum2:+f(actval+n*delta)
     1710       ENDIF
     1720     OTHERWISE
     1730       integral:=0
     1740     ENDCASE
     1750   ENDFOR
     1760   IF type$="S" THEN
     1770     RETURN (delta/6)*(f(lbound)+f(rbound)+4*sum1+2*sum2)
     1780   ELSE
     1790     RETURN integral
     1800   ENDIF
     1810 ENDFUNC
     1820 //
     1830 FUNC f1(x) CLOSED
     1840   RETURN x^3
     1850 ENDFUNC
     1860 //
     1870 FUNC f2(x) CLOSED
     1880   RETURN 1/x
     1890 ENDFUNC
     1900 //
     1910 FUNC f3(x) CLOSED
     1920   RETURN x
     1930 ENDFUNC
 