# This worksheet implements the integration package used in Chapter 15
# of Introduction to Scientific Programming by Joseph L. Zachary.
# 
# Create an empty package.
# 
> integration := table();
# 
# This animates how the rectangular method converges as it integrages
# integrand between lo and hi.  It creates nframes frames.  In the first
# there are two rectangles, and the number of rectangles doubles in each
# subsequent frame.
# 
> integration[animateRectangular] := proc (integrand, lo, hi, nframes)
>   local frame,
>         rect;
# 
# This creates a plot showing how func can be divided into rectangles. 
# Each rectangle is of width step, there are steps rectangles, and the
# range of interest goes from lo to hi.
# 
>  rect := proc (func, step, steps, lo, hi)
>    local n, x, onerect;
> 
>    onerect := (h, n, f, lo) -> ([lo+(n-1)*h, 0], [lo+(n-1)*h, f(lo +
> n*h)], [lo+n*h, f(lo + n*h)]);
> 
>    plot({func(x),
>          [seq(onerect(step, n, func, lo), n=1..steps), [lo+steps*step,
> 0]]},
>         x=lo..hi);
>  end;
> 
>  plots[display]([seq(rect(integrand, (hi-lo)/2^frame, 2^frame, lo,
> hi),
>                 frame=1..nframes)],
>                 insequence=true,
>                 title=`Rectangular Approximation`);
> end;
# 
# This animates how the trapezoidal method converges as it integrages
# integrand between lo and hi.  It creates nframes frames.  In the first
# there are two trapezoids, and the number of rectangles doubles in each
# subsequent frame.
# 
> integration[animateTrapezoidal] := proc (integrand, lo, hi, nframes)
>   local frame,
>         trap;
# 
# This creates a plot showing how func can be divided into trapezoids. 
# Each trapezoid is of width step, there are steps trapezoids, and the
# range of interest goes from lo to hi.
# 
>   trap := proc (func, step, steps, lo, hi)
> 
>     local n, x, onetrap;
> 
>     onetrap := (h, n, f, lo) -> ([lo+(n-1)*h, 0], [lo+(n-1)*h, f(lo +
> (n-1)*h)], [lo+n*h, f(lo + n*h)]);
> 
>     plot({func(x),
>           [seq(onetrap(step, n, func, lo), n=1..steps),
> [lo+steps*step, 0]]},
>          x=lo..hi);
> 
>     end;
> 
> 
>   plots[display]([seq(trap(integrand, (hi-lo)/2^frame, 2^frame, lo,
> hi),
>                  frame=1..nframes)],
>                  insequence=true,
>                  title= `Trapezoidal Approximation`);
> 
> end;
>