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.


// Rosetta Code problem: https://rosettacode.org/wiki/Numerical_integration
// by Jjuanhdez, 06/2022
 
print "function     range       steps  leftrect      midrect       rightrect     trap          simpson "
frmt$ = "%1.10f"
print "f(x) = x^3   0 - 1         100  ";
Integrate(0.0, 1.0, 1, 100)
print "f(x) = 1/x   1 - 100      1000  ";
Integrate(1.0, 100.0, 2, 1000)
frmt$ = "%8.3f"
print "f(x) = x     0 - 5000  5000000  ";
Integrate(0.0, 5000.0, 3, 5000000)
print "f(x) = x     0 - 6000  6000000  ";
Integrate(0.0, 6000.0, 3, 6000000)
end
 
sub Func(FN, X)         //Return F(X) for function number FN
    switch FN
    case 1
        return X ^ 3
    case 2
        return 1.0 / X
    case 3
        return X
    default
        return 0.0
    end switch
end sub
 
sub Integrate(A, B, FN, N)    //Display area under curve for function FN
    // A, B, FN    limits A, B, and number of slices N
 
    DX = (B-A)/N
    X = A 
    Area = 0.0             //rectangular left
    for i = 1 to N
        Area = Area + Func(FN,X)*DX 
        X = X + DX
    next i
    print str$(Area, frmt$);
    X = A
    Area = 0.0             //rectangular right
    for i = 1 to N
        X = X + DX 
        Area = Area + Func(FN,X)*DX
    next i
    print "  "; 
    print str$(Area, frmt$);
    X = A + DX / 2.0
    Area = 0.0             //rectangular mid point
    for i = 1 to N
        Area = Area + Func(FN,X)*DX 
        X = X + DX
    next i
    print "  "; 
    print str$(Area, frmt$);
    X = A 
    Area = 0.0             //trapezium
    for i = 1 to N
        Area = Area + (Func(FN,X)+Func(FN,X + DX))/2.0*DX 
        X = X + DX
    next i
    print "  "; 
    print str$(Area, frmt$);
    X = A 
    Area = 0.0             //Simpson's rule
    for i = 1 to N
        Area = Area + DX/6.0*(Func(FN,X) + 4.0*Func(FN,(X+X + DX)/2.0) + Func(FN,X + DX))
        X = X + DX
    next i
    print "  "; 
    print str$(Area, frmt$)
end sub