#The SUMS package. Created by Carmen Artino, January, 1993.
#This package consists of five PROCedures, LEFT, RIGHT, TRAP, SIMP and MID
#which compute numerical approximations to definite integrals.
#LEFT and RIGHT compute the left- and right-hand Riemann sums respectively.
#TRAP is the trapezoidal rule, SIMP is simpson's rule and MID is
#the midpoint rule.
#
#This package was designed to be used with the Harvard Calculus project and
#conforms to the notation used in that text; in particular, SIMP is the
#weighted average of MID and TRAP. This means that when using SIMP, n need
#not be even and SIMP(n) will give the same result as simpson(2*n).
#
#Note: This packages automatically loads in the student calculus package.
#
#Added in 1994: The PROCedures Left and Right which compute the symbolic
#left- and right-hand Riemann sums respectively.
#################################

`help/text/LEFT` := TEXT(
`FUNCTION: LEFT - computes the left hand Riemann sum of an expression `,
`defined over an interval.`,
``,
`CALLING SEQUENCE:`,
`  LEFT(expr,intrvl,n)`,
``,
`PARAMETERS:`,
`  expr   - an expression in a single variable.`,
`  intrvl - an interval of the form var = a..b where var is the variable `,
`           appearing in expr.`,
`  n      - a positive integer indicating the number of subintervals over `,
`           which to compute the left hand sum.`,
``,
`SYNOPSIS:`,
`  A typical call to LEFT will compute the Riemann sum of the`,
`  expression over the interval indicated using n subintervals`,
`  and present the result as a floating point number.`,
``,
``,
``,
`EXAMPLES:`,
`> LEFT(x^2, x=0..1, 4);`,
``,
`                         .2187500000`,
``,
`> LEFT(exp(-t^2), t=1..3, 5);`,
``,
`                         .2227872162`,
``,
`see also Left, RIGHT, Right`):

#################################

`help/text/Left` := TEXT(
`FUNCTION: Left - computes the symbolic left hand Riemann sum of`,
`an expression defined over an interval.`,
``,
`CALLING SEQUENCE:`,
`  Left(expr,intrvl,n)`,
``,
`PARAMETERS:`,
`  expr   - an expression in a single variable.`,
`  intrvl - an interval of the form var = a..b where var is the variable `,
`           appearing in expr.`,
`  n      - a positive integer indicating the number of subintervals over `,
`           which to compute the left hand sum.`,
``,
`SYNOPSIS:`,
`  A typical call to Left will compute the Riemann sum symbolically of the`,
`  expression over the interval indicated using n subintervals.`,
``,
``,
``,
`EXAMPLES:`,
`> Left(x^2, x=0..b, 4);`,
``,
`                                3`,
`                          7/32 b`,
``,
`simplify(Left(x^2, x=0..b, n));`,
``,
`                               3     2`,
`                              b  (2 n  - 3 n + 1)`,
`                          1/6 -------------------`,
`                                        2`,
`                                       n`,
``,
`limit(", n = infinity);`,
``,
`                                        3`,
`                                   1/3 b`,
``,
`Left(exp(-t^2), t=0..b, 3);`,
``,
`                                       2               2`,
`                 1/3 b (1 + exp(- 1/9 b ) + exp(- 4/9 b ))`,
``,
`see also LEFT, RIGHT, Right`):

#################################

`help/text/RIGHT` := TEXT(
`FUNCTION: RIGHT - computes the right hand Riemann sum of an expression `,
`defined over an interval.`,
``,
`CALLING SEQUENCE:`,
`  RIGHT(expr,intrvl,n)`,
``,
`PARAMETERS:`,
`  expr   - an expression in a single variable.`,
`  intrvl - an interval of the form var = a..b where var is the variable `,
`           appearing in expr.`,
`  n      - a positive integer indicating the number of subintervals over `,
`           which to compute the right hand sum.`,
``,
`SYNOPSIS:`,
`  A typical call to RIGHT will compute the Riemann sum of the`,
`  expression over the interval indicated using n subintervals`,
`  and present the result as a floating point number.`,
``,
``,
``,
`EXAMPLES:`,
`> RIGHT(x^2, x=0..1, 4);`,
``,
`                         .4687500000`,
``,
`> RIGHT(exp(-t^2), t=1..3, 5);`,
``,
`                         .07568480361`,
``,
``,
`see also Right, LEFT, Left`):

#################################

`help/text/Right` := TEXT(
`FUNCTION: Right - computes the symbolic right hand Riemann sum of`,
`an expression defined over an interval.`,
``,
`CALLING SEQUENCE:`,
`  Right(expr,intrvl,n)`,
``,
`PARAMETERS:`,
`  expr   - an expression in a single variable.`,
`  intrvl - an interval of the form var = a..b where var is the variable `,
`           appearing in expr.`,
`  n      - a positive integer indicating the number of subintervals over `,
`           which to compute the right hand sum.`,
``,
`SYNOPSIS:`,
`  A typical call to Right will compute the Riemann sum of the`,
`  expression over the interval indicated using n subintervals.`,
``,
``,
``,
`EXAMPLES:`,
`> Right(x^2, x=0..b, 4);`,
``,
`                                   15   3`,
`                                  ---- b`,
`                                   32`,
``,
`> Right(exp(-t^2), t=0..b, 3);`,
``,
`                               2               2           2`,
`             1/3 b (exp(- 1/9 b ) + exp(- 4/9 b ) + exp(- b ))`,
``,
``,
`see also RIGHT, LEFT, Left`):

#################################

`help/text/TRAP` := TEXT(
`FUNCTION: TRAP - numerically evaluates the definite integral of an `,
`expression defined over an interval.`,
``,
`CALLING SEQUENCE:`,
`  TRAP(expr,intrvl,n)`,
``,
`PARAMETERS:`,
`  expr   - an expression in a single variable.`,
`  intrvl - an interval of the form var = a..b where var is the variable `,
`           appearing in expr.`,
`  n      - a positive integer indicating the number of subintervals over `,
`           which to make the computation.`,
``,
`SYNOPSIS:`,
`  The function TRAP computes a numerical approximation to an integral`,
`  using the Trapezoidal rule.`,
``,
``,
``,
`EXAMPLES:`,
`> TRAP(x^2, x=0..1, 4);`,
``,
`                           .3437500000`,
``,
`> TRAP(exp(-t^2), t=1..3, 5);`,
``,
`                           .1492360099`,
``,
``,
`see also SIMP  MID`):

#################################

`help/text/SIMP` := TEXT(
`FUNCTION: SIMP - numerically evaluates the definite integral of an `,
`expression defined over an interval.`,
``,
`CALLING SEQUENCE:`,
`  SIMP(expr,intrvl,n)`,
``,
`PARAMETERS:`,
`  expr   - an expression in a single variable.`,
`  intrvl - an interval of the form var = a..b where var is the variable `,
`           appearing in expr.`,
`  n      - a positive integer indicating the number of subintervals `,
`           over which to make the computation. n need not be even`,
``,
`SYNOPSIS:`,
`  The function SIMP computes a numerical approximation to an integral`,
`  using a variation of Simpson's rule. SIMP is the weighted sum of TRAP`,
`  and MID`,
``,
``,
``,
`EXAMPLES:`,
`> SIMP(x^2, x=0..1, 4);`,
``,
`                          .3333333333`,
``,
`> SIMP(exp(-t^2), t=1..3, 6);`,
``,
`                          .1392862385`,
``,
``,
`see also TRAP  MID`):

#################################

`help/text/MID` := TEXT(
`FUNCTION: MID - numerical approximation to an integral`,
``,
`CALLING SEQUENCE:`,
`   MID(expr, intrvl, n)`,
``,
`PARAMETERS:`,
`   expr   - an algebraic expression in a single variable`,
`   intrvl - an interval of the form var = a..b where var is the variable `,
`            appearing in expr.`,
`   n      - a positive integer indicating the number of subintervals over`,
`            which to make the computation.`,
``,
`SYNOPSIS:`,
` The function MID computes a numerical approximation to a definite`,
` integral using rectangles.  The height of each rectangle (box) is`,
` determined by the value of the function at the midpoint of each`,
` interval.`,
``,
``,
``,
`EXAMPLES:`,
`> MID(x^2, x=0..1, 4);`,
``,
`                            .3281250000`,
``,
`> MID(exp(-t^2), t=1..3, 5);`,
``,
`                            .1344371675`,
``,
``,
`see also TRAP  SIMP`):

#################################

with(student):

#################################

LEFT := proc(expr, intrvl, n)
        local a, b, h, s, var, i;
        var := op(1,intrvl):
        a:=op(1,op(2,intrvl)): b:=op(2,op(2,intrvl)):
        h := (b-a)/n:
        s := 0:
        for i from 0 to n-1 do
          s := s + evalf(h*subs(var=a+i*h,expr)):
        od:
        s;
        end:

#################################

RIGHT := proc(expr, intrvl, n)
         local a, b, h, s, var, i;
         var := op(1,intrvl):
         a:=op(1,op(2,intrvl)): b:=op(2,op(2,intrvl)):
         h := (b-a)/n:
         s := 0:
         for i from 1 to n do
           s := s + evalf(h*subs(var=a+i*h,expr)):
         od:
         s;
         end:

#################################

MID :=   proc(expr, intrvl, n)
         local a, b, h, s, var, i;
         var := op(1,intrvl):
         a:=op(1,op(2,intrvl)): b:=op(2,op(2,intrvl)):
         h := (b-a)/n:
         s := 0:
         for i from 1 to n do
           s := s + evalf(h*subs(var=a+((2*i-1)/2)*h,expr)):
         od:
         s;
         end:

#################################

TRAP := proc(expr, intrvl, n)
        local s;
        s := (LEFT(expr, intrvl, n) + RIGHT(expr, intrvl, n))/2:
        s;
        end:

#################################

SIMP := proc(expr, intrvl, n)
         evalf((2*MID(expr,intrvl,n)+TRAP(expr,intrvl,n))/3);
        end:

#################################

Left := proc(expr, intrvl, n)
          value(leftsum(expr, intrvl, n));
        end:

#################################

Right := proc(expr, intrvl, n)
          value(rightsum(expr, intrvl, n));
        end:

#################################

lrbox := proc(expr, intrvl, n)
         local l, r, a, b, var;
         var := op(1,intrvl):
         a:=op(1,op(2,intrvl)): b:=op(2,op(2,intrvl)):
         l := leftbox(expr, var = a..b, n, color = black):
         r := rightbox(expr, var = a..b, n, color = red):
         plot({l,r});
end:

#################################