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.


Prototype.d TestFunction(Arg.d)
 
Procedure.d LeftIntegral(Start, Stop, Steps, *func.TestFunction)
  Protected.d n=(Stop-Start)/Steps, sum, x=Start
  While x <= Stop-n
    sum + n * *func(x)
    x + n
  Wend
  ProcedureReturn sum
EndProcedure
 
Procedure.d MidIntegral(Start, Stop, Steps, *func.TestFunction)
  Protected.d n=(Stop-Start)/Steps, sum, x=Start
  While x <= Stop-n
    sum + n * *func(x+n/2)
    x + n
  Wend
  ProcedureReturn sum
EndProcedure
 
Procedure.d RightIntegral(Start, Stop, Steps, *func.TestFunction)
  Protected.d n=(Stop-Start)/Steps, sum, x=Start
  While x < Stop
    x + n
    sum + n * *func(x)
  Wend
  ProcedureReturn sum
EndProcedure
 
Procedure.d Trapezium(Start, Stop, Steps, *func.TestFunction)
  Protected.d n=(Stop-Start)/Steps, sum, x=Start
  While x<=Stop
    sum + n * (*func(x) + *func(x+n))/2
    x+n
  Wend
  ProcedureReturn sum
EndProcedure
 
Procedure.d Simpson(Start, Stop, Steps, *func.TestFunction)
  Protected.d n=(Stop-Start)/Steps, sum1, sum2, x=Start
  Protected i
  For i=0 To steps-1
    sum1+ *func(Start+n*i+n/2)
  Next
  For i=1 To Steps-1
    sum2+ *func(Start+n*i)
  Next
  ProcedureReturn n * (*func(Start)+ *func(Stop)+4*sum1+2*sum2) / 6
EndProcedure
 
;- Set up functions to integrate
Procedure.d Test1(n.d)
  ProcedureReturn n*n*n  
EndProcedure
 
Procedure.d Test2(n.d)
  ProcedureReturn 1/n  
EndProcedure
 
; This function should be integrated as a integer function, but for 
; comparably this will stay as a float.
Procedure.d Test3(n.d)
  ProcedureReturn n
EndProcedure
 
;- Test the code & present the results
CompilerIf #PB_Compiler_Debugger
  MessageRequester("Notice!","Running this program in Debug-mode will be slow")
CompilerEndIf
 
; = 0.25
Define Answer$
Answer$="Left     ="+StrD(LeftIntegral (0,1,100,@Test1()))+#CRLF$
Answer$+"Mid      ="+StrD(MidIntegral  (0,1,100,@Test1()))+#CRLF$
Answer$+"Right    ="+StrD(RightIntegral(0,1,100,@Test1()))+#CRLF$
Answer$+"Trapezium="+StrD(Trapezium    (0,1,100,@Test1()))+#CRLF$
Answer$+"Simpson  ="+StrD(Simpson      (0,1,100,@Test1()))
MessageRequester("Answer should be 1/4",Answer$)
 
; = Ln(100) e.g. ~4.60517019...
Answer$="Left     ="+StrD(LeftIntegral  (1,100,1000,@Test2()))+#CRLF$
Answer$+"Mid      ="+StrD(MidIntegral   (1,100,1000,@Test2()))+#CRLF$
Answer$+"Right    ="+StrD(RightIntegral (1,100,1000,@Test2()))+#CRLF$
Answer$+"Trapezium="+StrD(Trapezium     (1,100,1000,@Test2()))+#CRLF$
Answer$+"Simpson  ="+StrD(Simpson       (1,100,1000,@Test2()))
MessageRequester("Answer should be Ln(100), e.g. ~4.60517019",Answer$)
 
; 12,500,000
Answer$="Left     ="+StrD(LeftIntegral  (0,5000,5000000,@Test3()))+#CRLF$
Answer$+"Mid      ="+StrD(MidIntegral   (0,5000,5000000,@Test3()))+#CRLF$
Answer$+"Right    ="+StrD(RightIntegral (0,5000,5000000,@Test3()))+#CRLF$
Answer$+"Trapezium="+StrD(Trapezium     (0,5000,5000000,@Test3()))+#CRLF$
Answer$+"Simpson  ="+StrD(Simpson       (0,5000,5000000,@Test3()))
MessageRequester("Answer should be 12,500,000",Answer$)
 
; 18,000,000
Answer$="Left     ="+StrD(LeftIntegral  (0,6000,6000000,@Test3()))+#CRLF$
Answer$+"Mid      ="+StrD(MidIntegral   (0,6000,6000000,@Test3()))+#CRLF$
Answer$+"Right    ="+StrD(RightIntegral (0,6000,6000000,@Test3()))+#CRLF$
Answer$+"Trapezium="+StrD(Trapezium     (0,6000,6000000,@Test3()))+#CRLF$
Answer$+"Simpson  ="+StrD(Simpson       (0,6000,6000000,@Test3()))
MessageRequester("Answer should be 18,000,000",Answer$) 