Lab 7: Student Assignment

Week 7, February 22

MAT 229, Spring 2021

Numerical Integration Examples

This lab begins by asking you to evaluate some “initialization code” provided by Prof. Al Hibbard.
Under the “Evaluation Tab” above you’ll see “Evaluate Notebook”: select that, and allow initialization. When you’re done, come back up here to the beginning.

Numerical Integration Code: thanks to Al Hibbard  -- HibbardA@central.edu

Prelims

(* Needs["BarCharts`"];
Needs["Histograms`"];
Needs["PieCharts`"];*)
Needs["Utilities`FilterOptions`"];

Off[General::spell,General::spell1];

Usage statements

ShowSeqSer::usage="ShowSeqSer[seq,{start,end}] uses
ListPlot to display both the underlying sequence and the
sequence of partial sums for series determined by seq. By
default, this shows every term from start to end. One can
use ShowSeqSer[seq,{start,end, step}] to use other steps.";

ShowArcLength::usage="ShowArcLength[fn, {var, min, max},
pts, opts] creates a plot of fn over the domain [min, max]
and then also draws a linear approximation using line
segments using x coordinates at values of pts. The length of
the line segment path is also given. The final parameter
opts can take any plot options. Note that one can use
ShowArcLength[fn, {var, min, max}, n, opts] and then n
points are uniformly chosen in the interval [min, max].";

circumscribed::usage="circumscribed[fn, {var, min,max},n,opts]
approximates the area under the curve given by fn and above the x-axis
over the domain [min,max] when n circumscribed rectangles are used.
circumscribed[{fn1, fn2}, {var, min,max},n,opts] acts similarly except
that the rectangles approximate the area between the curves fn1 and fn2.
The parameter n defaults to 8 if not specified and opts is for options.
One can add the option PtsInfoToo -> True (defaults to False) and the
pair {xpoints, area} will be returned where xpoints are the points in
the interval [min,max] that were used and area is the area of the
rectangles. Note this function returns the area in exact form; use
ncircumscribed for decimal form. Note also that the word circumscribed
is also a type used in ShowApproximations.";

circumscribedR::usage="circumscribedR[fn, {var, min,max},n,color]
creates the graphics to display the rectangular approximation of the
area under the curve fn (and above the x-axis) over the domain [min,max]
when n rectangles are used in the approximation. The color has a default
choice, so it can be omitted. The function returns a list {rlist, llist}
of graphics primitives necessary to draw the rectangles and their
perimeters (the lines list). circumscribedR[{fn1, fn2}, {var,
min,max},n,color] works similarly except that the rectangles are
designed to approximate the area between the curves fn1 and fn2.";

closest::usage="closest[{f1, f2}, {var1, a, b}, {var2, c, d}, {n1,
n2}, opts] returns the (exact) closest approximation (using the point
closest to the origin) for the volume between the functions of two
variables f1 and f2, over the rectangle [a, b] x [c, d] using n1*n2
divisions. Alternatively, closest[{f1, f2}, {var1, a, b}, {var2, c, d},
n, opts] assumes n1 = n2 = n. Also, closest[f, ...] can be used where f1
= 0 and f2 = f. The option PtsInfoToo (taking True or False, the default
being False) can be added also.";

closestC::usage="";

Color::usage="Color is an option for some functions, including
ShowApproximations. It is used to specify the color to be used in
showing an approximation. Color -> Red will color them Red.";

DisplayArea::usage="DisplayArea is an an option for the
ShowApproximations function which will cause the value of the area
estimate to be displayed if set to True (which it is by default).
DisplayArea -> False will suppress printing any area information.";

inscribed::usage="inscribed[fn, {var, min,max},n,opts] approximates
the area under the curve given by fn and above the x-axis over the
domain [min,max] when n inscribed rectangles are used. inscribed[{fn1,
fn2}, {var, min,max},n,opts] acts similarly except that the rectangles
approximate the area between the curves fn1 and fn2. The parameter n
defaults to 8 if not specified and opts is for options. One can add the
option PtsInfoToo -> True (defaults to False) and the pair {xpoints,
area} will be returned where xpoints are the points in the interval
[min,max] that were used and area is the area of the rectangles. Note
this function returns the area in exact form; use ninscribed for decimal
form. Note also that the word inscribed is also a type used in
ShowApproximations.";

inscribedR::usage="inscribedR[fn, {var, min,max},n,color] creates the
graphics to display the rectangular approximation of the area under the
curve fn (and above the x-axis) over the domain [min,max] when n
rectangles are used in the approximation. The color has a default
choice, so it can be omitted. The function returns a list {rlist, llist}
of graphics primitives necessary to draw the rectangles and their
perimeters (the lines list). inscribedR[{fn1, fn2}, {var,
min,max},n,color] works similarly except that the rectangles are
designed to approximate the area between the curves fn1 and fn2.";

left::usage="left[fn, {var, min,max},n,opts] approximates the area
under the curve given by fn and above the x-axis over the domain
[min,max] when n rectangles using left endpoints are used. left[{fn1,
fn2}, {var, min,max},n,opts] acts similarly except that the rectangles
approximate the area between the curves fn1 and fn2. The parameter n
defaults to 8 if not specified and opts is for options. One can add the
option PtsInfoToo -> True (defaults to False) and the pair {xpoints,
area} will be returned where xpoints are the points in the interval
[min,max] that were used and area is the area of the rectangles. Note
this function returns the area in exact form; use nleft for decimal
form. Note also that the word left is also a type used in
ShowApproximations.";

LeftApprox::usage="LeftApprox[fn,{var,min,max},n,opts] returns an
exact value for approximating the area under the function fn over the
interval [min, max] when n (which defaults to 8 if unspecified)
subintervals are used with the left-endpoint rectangle method. Adding
the option PtsInfoToo -> True (which defaults to False) will
additionally return the coordinates of the n points in the interval
[min, max] which are sampled. LeftApprox[{fn1,fn2},{var,min,max},n,opts]
is similar except that the area in question is between the two functions
fn1 and fn2.";

leftR::usage="leftR[fn, {var, min,max},n,color] creates the graphics
to display the rectangular approximation of the area under the curve fn
(and above the x-axis) over the domain [min,max] when n rectangles are
used in the approximation. The color has a default choice, so it can be
omitted. The function returns a list {rlist, llist} of graphics
primitives necessary to draw the rectangles and their perimeters (the
lines list). leftR[{fn1, fn2}, {var, min,max},n,color] works similarly
except that the rectangles are designed to approximate the area between
the curves fn1 and fn2.";

midpoint::usage="midpoint[fn, {var, min,max},n,opts] approximates the
area under the curve given by fn and above the x-axis over the domain
[min,max] when n rectangles using midpoints are used. midpoint[{fn1,
fn2}, {var, min,max},n,opts] acts similarly except that the rectangles
approximate the area between the curves fn1 and fn2. The parameter n
defaults to 8 if not specified and opts is for options. One can add the
option PtsInfoToo -> True (defaults to False) and the pair {xpoints,
area} will be returned where xpoints are the points in the interval
[min,max] that were used and area is the area of the rectangles. Note
this function returns the area in exact form; use nmidpoint for decimal
form. Note also that the word midpoint is also a type used in
ShowApproximations.\n\n midpoint[{f1, f2}, {var1, a, b}, {var2, c, d},
{n1, n2}, opts] returns the \ (exact) midpoint approximation for the
volume between the functions of two \ variables f1 and f2, over the
rectangle [a, b] x [c, d] using n1*n2 \ divisions. Alternatively,
midpoint[{f1, f2}, {var1, a, b}, {var2, c, d}, n, \ opts] assumes n1 =
n2 = n. Also, midpoint[f, ...] can be used where f1 = 0 \ and f2 = f.
The option PtsInfoToo (taking True or False, the default being \ False)
can be added also.\n\nmidpoint[{f1, f2}, {var1, a, b}, {var2, c, d}, {n1,
n2}, opts] returns the (exact) midpoint approximation for the volume
between the functions of two variables f1 and f2, over the rectangle [a,
b] x [c, d] using n1*n2 divisions. Alternatively, midpoint[{f1, f2},
{var1, a, b}, {var2, c, d}, n, opts] assumes n1 = n2 = n. Also,
midpoint[f, ...] can be used where f1 = 0 and f2 = f. The option
PtsInfoToo (taking True or False, the default being False) can be added
also.";

MidpointApprox::usage="MidpointApprox[fn,{var,min,max},n,opts] returns
an exact value for approximating the area under the function fn over the
interval [min, max] when n (which defaults to 8 if unspecified)
subintervals are used with the midpoint rectangle method. Adding the
option PtsInfoToo -> True (which defaults to False) will additionally
return the coordinates of the n points in the interval [min, max] which
are sampled. MidpointApprox[{fn1,fn2},{var,min,max},n,opts] is similar
except that the area in question is between the two functions fn1 and
fn2.";

midpointC::usge="";

midpointR::usage="midpointR[fn, {var, min,max},n,color] creates the
graphics to display the rectangular approximation of the area under the
curve fn (and above the x-axis) over the domain [min,max] when n
rectangles are used in the approximation. The color has a default
choice, so it can be omitted. The function returns a list {rlist, llist}
of graphics primitives necessary to draw the rectangles and their
perimeters (the lines list). midpointR[{fn1, fn2}, {var,
min,max},n,color] works similarly except that the rectangles are
designed to approximate the area between the curves fn1 and fn2.";

ncircumscribed::usage="ncircumscribed[fn, {var, min,max},n,opts] calls
on circumscribed[fn, {var, min,max},n,opts] and returns the value in
decimal form. Similarly for ncircumscribed[{fn1, fn2}, {var,
min,max},n,opts].";

nclosest::usage="nclosest[{f1, f2}, {var1, a, b}, {var2, c, d}, {n1,
n2}, opts] is similar to closest except that the returned value is given
in decimal form.";

ninscribed::usage="ninscribed[fn, {var, min,max},n,opts] calls on
inscribed[fn, {var, min,max},n,opts] and returns the value in decimal
form. Similarly for ninscribed[{fn1, fn2}, {var, min,max},n,opts].";

nleft::usage="nleft[fn, {var, min,max},n,opts] calls on left[fn, {var,
min,max},n,opts] and returns the value in decimal form. Similarly for
nleft[{fn1, fn2}, {var, min,max},n,opts].";

NLeftApprox::usage="NLeftApprox[fn,{var,min,max},n,opts] returns a
decimal value for approximating the area under the function fn over the
interval [min, max] when n (which defaults to 8 if unspecified)
subintervals are used with the left-endpoint rectangle method. Adding
the option PtsInfoToo -> True (which defaults to False) will
additionally return the coordinates of the n points in the interval
[min, max] which are sampled. LeftApprox[{fn1,fn2},{var,min,max},n,opts]
is similar except that the area in question is between the two functions
fn1 and fn2.";

nmidpoint::usage="nmidpoint[fn, {var, min,max},n,opts] calls on
midpoint[fn, {var, min,max},n,opts] and returns the value in decimal
form. Similarly for nmidpoint[{fn1, fn2}, {var, min,max},n,opts].\n\n
nmidpoint[{f1, f2}, {var1, a, b},
{var2, c, d}, {n1, n2}, opts] is similar to midpoint except that the
returned value is given in decimal form.";

NMidpointApprox::usage="NMidpointApprox[fn,{var,min,max},n,opts]
returns a decimal value for approximating the area under the function fn
over the interval [min, max] when n (which defaults to 8 if unspecified)
subintervals are used with the midpoint rectangle method. Adding the
option PtsInfoToo -> True (which defaults to False) will additionally
return the coordinates of the n points in the interval [min, max] which
are sampled. NMidpointApprox[{fn1,fn2},{var,min,max},n,opts] is similar
except that the area in question is between the two functions fn1 and
fn2.";

nrandom::usage="nrandom[fn, {var, min,max},n,opts] calls on random[fn,
{var, min,max},n,opts] and returns the value in decimal form. Similarly
for nrandom[{fn1, fn2}, {var, min,max},n,opts].";

nright::usage="nright[fn, {var, min,max},n,opts] calls on right[fn,
{var, min,max},n,opts] and returns the value in decimal form. Similarly
for nright[{fn1, fn2}, {var, min,max},n,opts].";

NRightApprox::usage="NRightApprox[fn,{var,min,max},n,opts] returns a
decimal value for approximating the area under the function fn over the
interval [min, max] when n (which defaults to 8 if unspecified)
subintervals are used with the right-endpoint rectangle method. Adding
the option PtsInfoToo -> True (which defaults to False) will
additionally return the coordinates of the n points in the interval
[min, max] which are sampled.
NRightApprox[{fn1,fn2},{var,min,max},n,opts] is similar except that the
area in question is between the two functions fn1 and fn2.";

nsimpson::usage="nsimpson[fn,{var,min,max},n,opts] returns a decimal
value for approximating the area under the function fn over the interval
[min, max] when n (which must be even) subintervals are used with
Simpson method. Adding the option PtsInfoToo -> True (which defaults to
False) will additionally return the coordinates of the n points in the
interval [min, max] which are sampled.";

NSimpsonApprox::usage="NSimpsonApprox[fn,{var,min,max},n] returns a
decimal value for approximating the area under the function fn over the
interval [min, max] when n (which defaults to 8 if unspecified)
subintervals are used with the Simpson method.";

ntrapezoid::usage="ntrapezoid[fn,{var,min,max},n,opts] returns a
decimal value for approximating the area under the function fn over the
interval [min, max] when n (which defaults to 8 if unspecified)
subintervals are used with trapezoid method. Adding the option
PtsInfoToo -> True (which defaults to False) will additionally return
the coordinates of the n points in the interval [min, max] which are
sampled.";

NTrapezoidApprox::usage="NTrapezoidApprox[fn,{var,min,max},n] returns
a decimal value for approximating the area under the function fn over
the interval [min, max] when n (which defaults to 8 if unspecified)
subintervals are used with the trapezoid method.
NTrapezoidApprox[{fn1,fn2},{var,min,max},n,opts] is similar except that
the area in question is between the two functions fn1 and fn2.";

PtsInfoToo::usage="PtsInfoToo is an option for several functions which
are used to determine approximations of area based on rectangles. If
PtsInfoToo -> True is added, the actual values of the points which in
the subinterval which were used to determine the heights will be also
given. Functions which can call this are left, right, midpoint, random,
inscribed, cirsumscribed and nxxxx variations of them.";

random::usage="random[fn, {var, min,max},n,opts] approximates the area
under the curve given by fn and above the x-axis over the domain
[min,max] when n rectangles using random points in the subinterval are
used. random[{fn1, fn2}, {var, min,max},n,opts] acts similarly except
that the rectangles approximate the area between the curves fn1 and fn2.
The parameter n defaults to 8 if not specified and opts is for options.
One can add the option PtsInfoToo -> True (defaults to False) and the
pair {xpoints, area} will be returned where xpoints are the points in
the interval [min,max] that were used and area is the area of the
rectangles. Note this function returns the area in exact form; use
nrandom for decimal form. Note also that the word random is also a type
used in ShowApproximations.";

randomR::usage="randomR[fn, {var, min,max},n,color] creates the
graphics to display the rectangular approximation of the area under the
curve fn (and above the x-axis) over the domain [min,max] when n
rectangles are used in the approximation. The color has a default
choice, so it can be omitted. The function returns a list {area,{rlist,
llist}} containing the estimated area and a list of graphics primitives
necessary to draw the rectangles and their perimeters (the lines list).
randomR[{fn1, fn2}, {var, min,max},n,color] works similarly except that
the rectangles are designed to approximate the area between the curves
fn1 and fn2.";

right::usage="right[fn, {var, min,max},n,opts] approximates the area
under the curve given by fn and above the x-axis over the domain
[min,max] when n rectangles using right endpoints are used. right[{fn1,
fn2}, {var, min,max},n,opts] acts similarly except that the rectangles
approximate the area between the curves fn1 and fn2. The parameter n
defaults to 8 if not specified and opts is for options. One can add the
option PtsInfoToo -> True (defaults to False) and the pair {xpoints,
area} will be returned where xpoints are the points in the interval
[min,max] that were used and area is the area of the rectangles. Note
this function returns the area in exact form; use nright for decimal
form. Note also that the word right is also a type used in
ShowApproximations.";

RightApprox::usage="RightApprox[fn,{var,min,max},n,opts] returns an
exact value for approximating the area under the function fn over the
interval [min, max] when n (which defaults to 8 if unspecified)
subintervals are used with the right-endpoint rectangle method. Adding
the option PtsInfoToo -> True (which defaults to False) will
additionally return the coordinates of the n points in the interval
[min, max] which are sampled.
RightApprox[{fn1,fn2},{var,min,max},n,opts] is similar except that the
area in question is between the two functions fn1 and fn2.";

rightR::usage="rightR[fn, {var, min,max},n,color] creates the graphics
to display the rectangular approximation of the area under the curve fn
(and above the x-axis) over the domain [min,max] when n rectangles are
used in the approximation. The color has a default choice, so it can be
omitted. The function returns a list {rlist, llist} of graphics
primitives necessary to draw the rectangles and their perimeters (the
lines list). rightR[{fn1, fn2}, {var, min,max},n,color] works similarly
except that the rectangles are designed to approximate the area between
the curves fn1 and fn2.";

ShowApproximations2::usage="ShowApproximations2[fn,{var,min,max},
type,n,opts] is similar to ShowApproximations. All graphics are in
Yellow and the PlotRange is not set to All.";

ShowApproximations::usage="ShowApproximations[fn,{var,min,max},type,n,opts] graphs the function fn
on the domain [min,max] and then places n (defaulting to 8) rectangles
between this curve and the x-axis. The variable type can be either left,
right, midpoint, random, inscribed, or circumscribed (and specifies how
the rectangles' heights are determined). Options include Color and
DisplayArea; each type of rectangle has a default color associated with
it but the option Color -> mycolor can be added. By default, DisplayArea
-> True is set, which controls whether the area estimate is displayed or
not. ShowApproximations[{fn1,fn2},{var,min,max},type,n,opts] is similar
except that the rectangles are drawn being bounded above by fn2 and
below by fn1. The parameter opts can also pass on any Graphics options
to the Plot command.\n\n ShowApproximations[{f1, f2}, {var1, a, b},
{var2, c, d}, type, {n1, n2}, opts] returns a graphic showing the
approximation determined by type for the volume between the functions of
two variables f1 and f2, over the rectangle [a, b] x [c, d] using n1*n2
divisions. Alternatively, ShowApproximations[f, ...] assumes f1 = 0 and
f2 = f.\n\nShowApproximations[{f1, f2}, {var1, a, b},
{var2, c, d}, type, {n1, n2}, opts] returns a graphic showing the
approximation determined by type for the volume between the functions of
two variables f1 and f2, over the rectangle [a, b] x [c, d] using n1*n2
divisions. Alternatively, ShowApproximations[f, ...] assumes f1 = 0 and
f2 = f.";

ShowSimpsonApprox::usage="ShowSimpsonApprox[fn,{var,min,max},n,opts]
is essentially ShowApproximations[fn,{var,min,max}, simpson,n,opts]
without PlotRange -> All.";

ShowTrapezoidApprox::usage="ShowTrapezoidApprox[fn,{var,min,max},
n,opts] is essentially ShowApproximations[fn,{var,min,max},
trapezoid,n,opts] without PlotRange -> All.";

simpson::usage="simpson[fn,{var,min,max},n,opts] returns an exact
value for approximating the area under the function fn over the interval
[min, max] when n (which must be even) subintervals are used with
Simpson method. Adding the option PtsInfoToo -> True (which defaults to
False) will additionally return the coordinates of the n points in the
interval [min, max] which are sampled.";

SimpsonApprox::usage="SimpsonApprox[fn,{var,min,max},n,opts] returns
an exact value for approximating the area under the function fn over the
interval [min, max] when n (which defaults to 8 if unspecified)
subintervals are used with the Simpson method.";

simpsonR::usage="simpsonR[fn,{var,min,max},n,opts] returns the
graphics for approximating the area under the function fn over the
interval [min, max] when n (which must be even) subintervals are used
with the Simpson method. Show[simpsonR[fn, {var,min,max},n],
DisplayFunction -> $DisplayFunction] will show the graphics used.";

TraceParametricPlot::usage="TraceParametricPlot[{f, g}, {t, t0, t1,
n}, opts] creates a series of parametric plots of the parametric
functions x = f[t] and y = g[t], starting at t = t0 and ending at t = t1
using n snapshots (defaulting to 8 if not specified). This is ready for
animation after being created. Any graphics options can be used for
opts.";

trapezoid::usage="trapezoid[fn,{var,min,max},n,opts] returns an exact
value for approximating the area under the function fn over the interval
[min, max] when n (which defaults to 8 if unspecified) subintervals are
used with trapezoid method. Adding the option PtsInfoToo -> True (which
defaults to False) will additionally return the coordinates of the n
points in the interval [min, max] which are sampled.";

TrapezoidApprox::usage="TrapezoidApprox[fn,{var,min,max},n,opts]
returns an exact value for approximating the area under the function fn
over the interval [min, max] when n (which defaults to 8 if unspecified)
subintervals are used with the trapezoid method.
TrapezoidApprox[{fn1,fn2},{var,min,max},n,opts] is similar except that
the area in question is between the two functions fn1 and fn2.";

trapezoidR::usage="trapezoidR[fn, {var,min,max},n,color] returns the
graphics primitives for approximating the area under the function fn
over the interval [min, max] when n (which defaults to 8 if unspecified)
subintervals are used with the trapezoid method. The color defaults to
Green if unspecified. Show[Graphics[trapezoidR[fn,
{var,min,max},n,color]]] will show the graphics used.";

Visual::usage="Visual is a value for the Mode option. Some type of
visualization will be given, in addition to a computation.";

RainbowPolarPlot::usage="RainbowPolarPlot[r,{var,min,max}, n,
opts] plots the polar function given by r and colors the graph
in a rainbow fashion (using Hue) using n different colors, one
for each (max-min)/n section. This allows one to see the path that
the curve takes as var is increased from min to max. One can
add optional Plot options, specified in opts.";

TracePolarPlot::usage="TracePolarPlot[r,{var,min,max},n,opts]
creates a series of n plots showing the graph of the polar function
successively being built up from var = min to var = max. The radius
of the curve at each break point is also shown. One can
add optional Plot options, specified in opts.
TracePolarPlot[r,{var,min,max},pts,opts] is similar except the
breakpoints are specified in the list pts.";

AverageValue::usage="AverageValue[fn,{var,min,max},opts]
calls FilledPlot with the function fn over [min,max] and
then overlays it with a rectangle with base the x axis and
top at height corresponding to the average value. The
average value is given as output. The parameter opts can
handle any Graphics options.";

Both::usage="Both is a value for the Side option. When
Side -> Both is given in ShowSecantToTan, ShowLimit, or
TableLimit, then animation or table begins from first the
left and then from the right.";

CutCardboard::usage="CutCardboard[width, length, tocut] is
similar to ShowCardboardBox except that it also shows (along
side) a depiction of how the resulting box appears (if cuts
are made at tocut).";

FilledPlotWithSquares::usage="FilledPlotWithSquares[fn,{var,min,max},size,opts] calls
FilledPlot with the function fn plotted over the domain
[min,max] and then overlays it with gridlines composed of
squares whose width is given by size. The parameter opts can
handle any Graphics options.";

GridCircle::usage="GridCircle[rad,wdth,opts] will draw a
circle with radius rad on a background of gridlines forming
squares whose width is given by wdth. Any Graphics options
can be added in opts.";

GrowAreaFunction::usage="GrowAreaFunction[fn,{var,min,max},n,opts] is a Do loop
which calls ShowAreaFunction n times, equally spaced, in the
interval [min,max], showing how this area function is built.
The parameter opts can handle any Graphics options.";

Left::usage="Left is a value for the Side option. When
Side -> Left is given in ShowSecantToTan, ShowLimit, or
TableLimit, then animation or table begins from the left.";

LimitTable::usage="LimitTable[fn, {var, min,
max},limitpt,n,opts] creates a table of two columns. The
first column has n values of var from var = min to var =
limitpt and then n more from var = limitpt to var = max. The
second column consists of the expression fn evaluated at
each of these values. The points in the first column become
exponentially closer to limitpt. The parameter n defaults to
7 if not given a different value. Optionally in opts you can
specify Side -> Right, Side -> Left or (by default) Side ->
Both, to indicate from which side(s) we approach limitpt.";

MakeOneParabola::usage="MakeOneParabola[fn, {var, first,
second}, opts] returns the graphics showing the unique
parabola which goes through the points (first,fn(first)),
((first+second)/2, fn((first+second)/2)), and (second,
fn(second)). This is used for the Simpson-related
functions.";

MarkIntervals::usage="MarkIntervals[fn, {var,min,max},
options] plots fn, fn' and fn'' and then draws vertical
lines where fn' = 0 and where fn'' = 0. See
PlotWithTangents.";

mesg1:=Print[StringJoin["The Blue dots are points whose
x-coordinate is a zero for the first derivative."]];

mesg2:=Print[StringJoin["The Magenta dots are points whose
x-coordinate is a zero for the second derivative."]];

mesg3:=Print[StringJoin["The Blue dots are points whose
x-coordinate is a zero for the first derivative and the
Magenta dots correspond to zeros of the second
derivative."]];

Side::usage="Side is the name of an
option of the function 	ShowSecantToTan. Values can be
specified as Left or Right; any other value causes a
default algorithm to decide from which side to illustrate
the progression of lines.";

NewtonSlideShow::usage="";

PlotWithTangents::usage="PlotWithTangents[fn,
{var,min,max}, list,options] assumes that the function fn is
in the variable var, plots it over the domain [min,max] and
then places tangents at e ach x value in list. One can add
any Plot options at the end. Additionally, one can set the
option ShowSlopes to True or False (defaults to True)
depending on whether or not the slopes of the tangent
line(s) should be given in the output. PlotWithTangents[fn,
{var,min,max}, n,options] places n tangent lines spread
uniformly across the interval [min, max].";

RectToTrap::usage="RectToTrap shows a 'Proof Without
Words' of how the midpoint approximation can be thought of
as a trapezoidal approximation.";

Right::usage="Right is a value for the Side option. When
Side -> Right is given in ShowSecantToTan, then animation
begins from the right.";

ShowArcLength::usage="ShowArcLength[fn, {var, min, max},
pts, opts] creates a plot of fn over the domain [min, max]
and then also draws a linear approximation using line
segments using x coordinates at values of pts. The length of
the line segment path is also given. The final parameter
opts can take any plot options. Note that one can use
ShowArcLength[fn, {var, min, max}, n, opts] and then n
points are uniformly chosen in the interval [min, max].";

ShowAreaFunction::usage="ShowAreaFunction[fn,{var,min,max},n,opts] calls FilledPlot
and graphs the function fn from min to max. Additionally, it
plots points of the area function on top of the FilledPlot.
By default, 8 points will of the area function will be
determined (equally spaced in the domain) unless n has a
different value. Another version of this function by the
same name allows one to insert a list of x coordinates for
the area function. This parameter is in the same position as
the n. The parameter opts can handle any Graphics options.";

ShowBetween::usage="ShowBetween[{fn1,fn2},{var,min,max},type,n,opts] actually
calls
ShowApproximations[{fn1,fn2},{var,min,max},type,n,opts]. See
ShowApproximations for information.";

ShowCardboardBox::usage="ShowCardboardBox[width, length,
tocut] displays a rectangluar region with dimensions width
by length and dotted lines representing where to cut the box
tocut inside the bounding rectangle.";

ShowChurchWindow::usage="ShowChurchWindow[radius,perimeter] illustrates a window
with perimeter units around surmounted with semi-circular
portion with the given radius.";

ShowChurchWindowDetails::usage="ShowChurchWindowDetails[radius,perimeter] is similar to
ShowChurchWindow with additional deails and hints given.";

ShowCuts::usage="ShowCuts[fn, {var, min, max}, pts, opts]
creates a FilledPlot between fn and -fn on the domain [min,
max] and then also places vertical lines (simulating were
cuts might be made) at x coordinates at values of pts. The
final parameter opts can take any plot options. Note that
one can use ShowCuts[fn, {var, min, max}, n, opts] and then
n points are uniformly chosen in the interval [min, max].";

ShowEpsilonDelta::usage="ShowEpsilonDelta[fn,{var,min,max},
limitpoint,epsilon,opts] graphs the function fn over the
interval [min,max], centered at var = limitpoint. Given
epsilon, appropriate lines are drawn at L + epsilon and L -
epsilon where L is the limit of fn as var approaches
limitpoint. Lines corresponding to limitpoint - d and
limitpoint + d are also drawn, where d is the corresponding
(maximum) delta for this scenario. Any Plot options can be
added for opts.";

ShowFence::usage="ShowFence[width, fencing] illustrates
how the field would appear when the given width is used, out
of the total amount (given by fencing).";

ShowFunctionAndDerivative::usage="ShowFunctionAndDerivative[fn, {var,min,max},list,options]
first calls PlotWithTangents and then produces a joined
ListPlot graph below it constructed from the Out of
PlotWithTangents (ie., points on the derivative of fn). This
function can also have list be replaced with an integer. See
PlotWithTangents.";

ShowFunctionAndDerivativeTogether::usage="ShowFunctionAndDerivativeTogether[fn, {var,min,max}, list,
options] is similar to ShowFunctionAndDerivative (see it)
except that both plots are shown in one plot rather than
two. Also see PlotWithTangents.";

ShowGuessedArea::usage="ShowGuessedArea[fn,{var,min,max},h,opts] calls FilledPlot
with the function fn over [min,max] and then overlays it
with a rectangle with base the x axis and top at height h.
The area to this is intended as a guess for the area under
the curve. The guessed area is calculated and a message is
given whether it is too high or low. The parameter opts can
handle any Graphics options.";

ShowLeftRightDifference::usage="ShowLeftRightDifference[fn, {var,min,max},n] shows the left
and right approximations superimposed on each other. The
number of subdivisions defaults to 8 if not specified.";

ShowLimit::usage="ShowLt[fn, {var, min, max}, limitpoint,
n, opts] runs a Do loop which calls ShowSecantTan. The loop
runs n times, which by default is 6 if not given, and the
points become exponentially closer to limitpoint Any Plot
options can be added optionally in opts.";

ShowLissajous::usage="ShowLissajous[f,g,pts] plots the
parametric equations x = f(t), y = g(t) for t in [0, 2Pi]
with the points pts labeled according to color scheme
indicating the order in which they are drawn.";

ShowMinMax::usage="ShowMinMax[fn,{var,min,max},opts] plots
the function fn with FilledPlot over the domain [min,max]
and then draws a horizontal line at the global maximum and
global minimum of this function on the given domain. The
parameter opts can handle any Graphics options.";

ShowPartialSum::usage="ShowPartialSum[fn, {var, min, max}]
assumes fn is a function of var and is sequence used for
defining a series. This will show the partial sum in an
unevaluated state from term var = min to var = max.";

ShowRectangle::usage="ShowRectangle[width,radius]
illustrates a rectangle with the given width can be
inscribed in a semi-circle with the given radius.";

ShowSecant::usage="ShowSecant[fn, {var, min, max}, pt1,
pt2, opts] will plot the function fn in the variable var on
the domain [min, max]. Added will be the secant from pt1 to
pt2. Any Plot options can be added optionally in opts.";

ShowSecantTan::usage="ShowSecantTan[fn, {var, min, max},
pt1, pt2, opts] will plot the function fn in the variable
var on the domain [min, max]. The points pt1 and pt2 can
either be just the x-values or in the list form {x,y}. Added
will be the secant from pt1 to pt2 and the tangent line to
fn at pt1. Any Plot options can be added optionally in
opts.";

ShowSecantToTan::usage="ShowSecantTan[fn, {var, min, max},
limitpoint, otherpoint, opts] will plot the function fn in
the variable var on the domain [min, max]. Added will be the
secant from limitpoint to otherpoint and the tangent line to
fn at limitpoint Any Plot options can be added optionally in
opts.";

ShowSecantToTan::usage="ShowSecantTan[fn, {var, min, max},
pt1, n, opts] runs a Do loop which calls ShowSecantTan. The
loop runs n times, which by default is 7 if not given, and
the points become exponentially closer to pt1. Any Plot
options can be added optionally in opts.";

ShowSeqSer::usage="ShowSeqSer[seq,{start,end}] uses
ListPlot to display both the underlying sequence and the
sequence of partial sums for series determined by seq. By
default, this shows every term from start to end. One can
use ShowSeqSer[seq,{start,end, step}] to use other steps.";

ShowSeries::usage="ShowSeries[ser,n] shows a graphical
display of the accumulating sequence of partial sums using
Hue to discriminate between the contributions of successive
terms of the underlying sequence. By examining the width of
successive terms, one can visualize whether a series is
likely to converge or not. By default, this starts with the
first term. One can use ShowSeries[ser,n,start] to start at
another point.";

ShowSlope::usage="ShowSlope[f,g,pts] shows the slope of
the parametric plot x = f(t), y = g(t) for t in [0, 2Pi]
with the points pts labeled according to color scheme
indicating the order in which they are drawn.";

ShowSlopes::usage="ShowSlopes is an option for
PlotWithTangents. When set to True, it shows the slopes of
the drawn tangent line(s) as the output.";

ShowSpeed::usage="ShowSpeed[f,g,pts] shows the 'speed' of
the parametric plot x = f(t), y = g(t) for t in [0, 2Pi]
with the points pts labeled according to color scheme
indicating the order in which they are drawn.";

ShowTrapError::usage="ShowTrapError[fn,{var,min,max},n,opts] illustrates the
error in using the trapezoid method by coloring the missed
region(s).";

ShowTurns::usage="ShowTurns[fn,{var,min,max},opts] chooses
'good' places to show tangents and then calls
PlotWithTangents to do the work.";

Side::usage="Side is the name of an option of the
functions 	ShowSecantToTan, ShowLimit, and TableLimit.
Values can be 	specified as Left, Right or Both.";
Right::usage="Right is a value for the Side option. When
Side -> Right is given in ShowSecantToTan, ShowLimit, or
TableLimit, then animation or table begins from the right.";

SomeFigures::usage="SomeFigures[] simply draws a number of
figures using Graphics primitives.";

SSfunc::usage="SSfunc[var,n] assumes n is a nine-digit
number (such as a social security number) and returns a
polynomial in the variable var using the digits of n as
coefficients, interweaving minus signs.";

Tangraph::usage="Tangraph[fn, {var,min,max}, value,
options] plots the function fn, draws the tangent to fn at x
= value, and also draws the graph of the derivative on the
interval [min, value]. See PlotWithTangents.";

TraceParametricPlot::usage="TraceParametricPlot[{f,g},{t,t0,t1,dt},opts] 	displays a
sequence of plots of the parametric equations x = f(t),
	y = g(t) for t starting at t0 and up to t1, incrementing
in units 	of (t1-t0)/dt. Note that dt defaults to 8, so
it can be omitted. Any 	graphics options can be used for
opts.";

EulerApproxList::usage="EulerApproxList[f, {t,Y}, h,n]";

SlopeField::usage="
	SlopeField[f,{t,t0,t1,dt},{y,y0,y1,dy},opts]";

Approximations

dist[pt1_List,pt2_List]:=Sqrt[(pt1[[1]]-pt2[[1]])^2+(pt1[[2]]-pt2[[2]])^2]

ShowArcLength[fn_,{var_,min_,max_},n_Integer,opts___]:=ShowArcLength[fn,{var,min,max},Table[i,{i,min,max,(max-min)/n}],opts]

ShowArcLength[fn_,{var_,min_,max_},xpts_List,opts___]:=Module[{crds,lns,est,act,fnx,n},fnx=Function[var,fn];
crds=Transpose[{xpts,Map[fnx,xpts]}];
n=Length[crds];
pts=Map[Point,crds];
lns=Line[crds];
est=Sum[dist[crds[[i]],crds[[i+1]]],{i,1,n-1}]//N;
Plot[fn,{var,min,max},Epilog->{{Black,PointSize[0.025],pts},{Blue,lns}},PlotRange->All,PlotLabel->StringJoin["(",ToString[n-1],") estimate = ",ToString[N[est,5]]],opts]]

modmidpointR[fn_,{var_,min_,max_},n_:8,color_:Magenta]:=modmidpointR[{0,fn},{var,min,max},n,color]

modmidpointR[{fn1_,fn2_},{var_,min_,max_},n_:8,color_:Magenta]:=Module[{xpoints,dx,i,lhgts1,lhgts2,rhgts1,rhgts2,hgts1,hgts2,ollist,llist,rlist,mdpts,slopes1, slopes2,deriv1,deriv2},dx=(max-min)/n;
xpoints=Table[i,{i,min,max,dx}];
mdpts=(Drop[xpoints,{1}]+Drop[xpoints,{n+1}])/2;
deriv1 =D[fn1,var];
deriv2 = D[fn2,var];
slopes1 = Map[(deriv1 /.var-> #)&,mdpts];
slopes2 = Map[(deriv2 /.var-> #)&,mdpts];
lhgts1 =Table[slopes1[[k]]*(mdpts[[k]]-dx/2)-slopes1[[k]]*mdpts[[k]]+fn1/.var->mdpts[[k]],{k,1, Length[mdpts]}];
lhgts2 =Table[slopes2[[k]]*(mdpts[[k]]-dx/2)-slopes2[[k]]*mdpts[[k]]+fn2/.var->mdpts[[k]],{k,1, Length[mdpts]}];
rhgts1 =Table[slopes1[[k]]*(mdpts[[k]]+dx/2)-slopes1[[k]]*mdpts[[k]]+fn1/.var->mdpts[[k]],{k,1, Length[mdpts]}];
rhgts2 =Table[slopes2[[k]]*(mdpts[[k]]+dx/2)-slopes2[[k]]*mdpts[[k]]+fn2/.var->mdpts[[k]],{k,1, Length[mdpts]}];
hgts1=Map[Function[var,fn1],mdpts];
hgts2=Map[Function[var,fn2],mdpts];
ollist={Blue,Table[{Line[{{xpoints[[i]],hgts2[[i]]},{xpoints[[i+1]],hgts2[[i]]}}],Line[{{xpoints[[i+1]],hgts2[[i]]},{xpoints[[i+1]],hgts1[[i]]}}],Line[{{xpoints[[i+1]],hgts1[[i]]},{xpoints[[i]],hgts1[[i]]}}],Line[{{xpoints[[i]],hgts1[[i]]},{xpoints[[i]],hgts2[[i]]}}]},{i,1,n}]};
llist = Table[Line[{{xpoints[[k]],lhgts2[[k]]},{xpoints[[k]],lhgts1[[k]]},{xpoints[[k+1]],rhgts1[[k]]},{xpoints[[k+1]],rhgts2[[k]]},{xpoints[[k]],lhgts2[[k]]}}],{k,1,n}];
rlist = {color,Table[Polygon[{{xpoints[[k]],lhgts2[[k]]},{xpoints[[k]],lhgts1[[k]]},{xpoints[[k+1]],rhgts1[[k]]},{xpoints[[k+1]],rhgts2[[k]]},{xpoints[[k]],lhgts2[[k]]}}],{k,1,n}]};
{rlist, ollist,llist}
]

Options[ShowApproximations]={DisplayArea->True,Color->RGBColor[0,0,1]};

ShowApproximations[{fn1_,fn2_},{var_,min_,max_},type_,n_:8,opts___]:=Module[{k,rlist,llist,randlist,prt,thex,deltax,mycolors,fnx,xxx,gr1,gr2,area,grlist,daQ,rc,finopts={opts}},daQ=DisplayArea/.{opts}/.Options[ShowApproximations];
rc=Color/.{opts}/.Options[ShowApproximations];
finopts=DeleteCases[finopts,Color->_];
finopts=DeleteCases[finopts,DisplayArea->_];
AppendTo[finopts,Axes->True];
mycolors=Switch[type,left,ColorData["Legacy", "Burlywood"],right,Blue,midpoint,Magenta,modmidpoint, Magenta,random,Red,inscribed,Cyan,circumscribed,Yellow,trapezoid,Green,simpson,Pink];
If[rc=!=RGBColor[0,0,1],mycolors=rc];
{area,grlist}=Switch[type,left,{nleft[{fn1,fn2},{var,min,max},n],leftR[{fn1,fn2},{var,min,max},n,mycolors]},right,{nright[{fn1,fn2},{var,min,max},n],rightR[{fn1,fn2},{var,min,max},n,mycolors]},midpoint,{nmidpoint[{fn1,fn2},{var,min,max},n],midpointR[{fn1,fn2},{var,min,max},n,mycolors]},
modmidpoint,{nmidpoint[{fn1,fn2},{var,min,max},n],modmidpointR[{fn1,fn2},{var,min,max},n,mycolors]},random,randomR[{fn1,fn2},{var,min,max},n,mycolors],inscribed,{ninscribed[{fn1,fn2},{var,min,max},n],inscribedR[{fn1,fn2},{var,min,max},n,mycolors]},circumscribed,{ncircumscribed[{fn1,fn2},{var,min,max},n],circumscribedR[{fn1,fn2},{var,min,max},n,mycolors]},trapezoid,{ntrapezoid[fn2,{var,min,max},n],trapezoidR[fn2,{var,min,max},n,mycolors]},simpson,{nsimpson[fn2,{var,min,max},n],simpsonR[fn2,{var,min,max},n,mycolors]}];
gr1=If[type===simpson,grlist,Graphics[{Thickness[0.005],grlist},Axes->True]];
gr2=Plot[{fn1,fn2},{var,min,max},Evaluate[finopts/.{gh_}->gh],PlotStyle->{{Black,Thickness[0.01]}},PlotLabel->If[daQ,Style[StringJoin["(",ToString[n],") est. (net) area: ",ToString[N[area,6]]],Smaller],""],Axes->True,DisplayFunction->Identity,PlotRange->All];
Show[{gr2,gr1,gr2},finopts,DisplayFunction->$DisplayFunction]]
(* added finopts at end on 3/13/05 *)

RectToTrap:=Show[Graphics[{{Yellow,Rectangle[{0,0},{1,1}]},{Red,Rectangle[{0,1},{1,1.2`}]},{Yellow,Rectangle[{1.5`,0},{2.5`,1}]},{Red,Polygon[{{1.5`,1},{2.5`,1},{2.5`,1.4`},{1.5`,1}}]},{Black,PointSize[0.01`],Point[{0.5`,1.2`}],Point[{2,1.2`}],Line[{{0,1},{1,1.4`}}],Line[{{1.5`,1.2`},{2.5`,1.2`}}]}}]]

ShowTrapError[fn_,{var_,min_,max_},n_:8,opts___]:=Module[{k,rlist,llist,prt,thex,deltax,mycolors,fnx,xxx,gr1,gr2},gr1=Plot[fn,{var,min,max},opts,DisplayFunction->Identity,Filling->Axis, FillingStyle->Red];gr2=ShowApproximations[fn,{var,min,max},trapezoid,n,opts,DisplayFunction->Identity];Show[{gr1,gr2},opts,DisplayFunction->$DisplayFunction]]

circumscribedR[fn_,{var_,min_,max_},n_:8,color_:Cyan]:=circumscribedR[{0,fn},{var,min,max},n,color]
circumscribedR[{fn1_,fn2_},{var_,min_,max_},n_:8,color_:Cyan]:=Module[{points1,points2,xpoints,area,dx,i,hgts1,hgts2,hgts,llist,rlist,rpts},dx=((max-min)/n);{points1,points2,area}=ncircumscribed[{fn1,fn2},{var,min,max},n,PtsInfoToo->True];xpoints=N[Table[i,{i,min,max,dx}]];hgts1=Last[Transpose[points1]];hgts2=Last[Transpose[points2]];hgts=Transpose[{hgts1,hgts2}];hgts=Transpose[(If[#1[[1]]>#1[[2]],{#1[[2]],#1[[2]]},#1]&)/@hgts];{hgts1,hgts2}=hgts;llist=Table[{Line[{{xpoints[[i]],hgts2[[i]]},{xpoints[[i+1]],hgts2[[i]]}}],Line[{{xpoints[[i+1]],hgts2[[i]]},{xpoints[[i+1]],hgts1[[i]]}}],Line[{{xpoints[[i+1]],hgts1[[i]]},{xpoints[[i]],hgts1[[i]]}}],Line[{{xpoints[[i]],hgts1[[i]]},{xpoints[[i]],hgts2[[i]]}}]},{i,1,n}];rlist={color,Table[Rectangle[{xpoints[[i]],hgts2[[i]]},{xpoints[[i+1]],hgts1[[i]]}],{i,1,n}]};{rlist,llist}]

circumscribed[fn_,{var_,min_,max_},n_:8,opts___]:=circumscribed[{0,fn},{var,min,max},n,opts]

circumscribed[{fn1_,fn2_},{var_,min_,max_},n_:8,opts___]:=Module[{points1,points2,dx,i,hgts,area,ptstoo},ptstoo=PtsInfoToo/.{opts}/.Options[circumscribed];dx=((max-min)/n);points1=If[fn1===0,Table[{i,0},{i,min+dx,max,dx}],If[PolynomialQ[fn1,var],Table[First[PolyMinMax[fn1,{var,min+(i-1) dx,min+i dx}]],{i,1,n}],Table[First[OtherMinMax[fn1,{var,min+(i-1) dx,min+i dx}]],{i,1,n}]]];points2=If[fn2===0,Table[{i,0},{i,min+dx,max,dx}],If[PolynomialQ[fn2,var],Table[Last[PolyMinMax[fn2,{var,min+(i-1) dx,min+i dx}]],{i,1,n}],Table[Last[OtherMinMax[fn2,{var,min+(i-1) dx,min+i dx}]],{i,1,n}]]];hgts1=Transpose[points1][[2]];hgts2=Transpose[points2][[2]];hgts=hgts2-hgts1/.x_?Negative->0;area=Plus@@(hgts dx);If[ptstoo,{points1,points2,area},area]]

inscribedR[fn_,{var_,min_,max_},n_:8,color_:Cyan]:=inscribedR[{0,fn},{var,min,max},n,color]

inscribedR[{fn1_,fn2_},{var_,min_,max_},n_:8,color_:Cyan]:=Module[{points1,points2,xpoints,area,dx,i,hgts1,hgts2,hgts,llist,rlist,rpts},dx=((max-min)/n);{points1,points2,area}=ninscribed[{fn1,fn2},{var,min,max},n,PtsInfoToo->True];xpoints=N[Table[i,{i,min,max,dx}]];hgts1=Last[Transpose[points1]];hgts2=Last[Transpose[points2]];hgts=Transpose[{hgts1,hgts2}];hgts=Transpose[(If[#1[[1]]>#1[[2]],{#1[[2]],#1[[2]]},#1]&)/@hgts];{hgts1,hgts2}=hgts;llist=Table[{Line[{{xpoints[[i]],hgts2[[i]]},{xpoints[[i+1]],hgts2[[i]]}}],Line[{{xpoints[[i+1]],hgts2[[i]]},{xpoints[[i+1]],hgts1[[i]]}}],Line[{{xpoints[[i+1]],hgts1[[i]]},{xpoints[[i]],hgts1[[i]]}}],Line[{{xpoints[[i]],hgts1[[i]]},{xpoints[[i]],hgts2[[i]]}}]},{i,1,n}];rlist={color,Table[Rectangle[{xpoints[[i]],hgts2[[i]]},{xpoints[[i+1]],hgts1[[i]]}],{i,1,n}]};{rlist,llist}]

inscribed[fn_,{var_,min_,max_},n_:8,opts___]:=inscribed[{0,fn},{var,min,max},n,opts]

inscribed[{fn1_,fn2_},{var_,min_,max_},n_:8,opts___]:=Module[{points1,points2,dx,i,hgts,area,ptstoo},ptstoo=PtsInfoToo/.{opts}/.Options[inscribed];dx=((max-min)/n);points1=If[fn1===0,Table[{i,0},{i,min+dx,max,dx}],If[PolynomialQ[fn1,var],Table[Last[PolyMinMax[fn1,{var,min+(i-1) dx,min+i dx}]],{i,1,n}],Table[Last[OtherMinMax[fn1,{var,min+(i-1) dx,min+i dx}]],{i,1,n}]]];points2=If[fn2===0,Table[{i,0},{i,min+dx,max,dx}],If[PolynomialQ[fn2,var],Table[First[PolyMinMax[fn2,{var,min+(i-1) dx,min+i dx}]],{i,1,n}],Table[First[OtherMinMax[fn2,{var,min+(i-1) dx,min+i dx}]],{i,1,n}]]];hgts1=Transpose[points1][[2]];hgts2=Transpose[points2][[2]];hgts=hgts2-hgts1/.x_?Negative->0;area=Plus@@(hgts dx);If[ptstoo,{points1,points2,area},area]]

OtherMinMax[fn_,{var_,min_,max_}]:=Module[{mins,maxs,d,both,fnx,hghts,nmin=N[min],nmax=N[max]},Off[FindMinimum::"fmlim",FindMinimum::"regex",FindMinimum::"fmns",Transpose::"nmtx",FindMinimum::"fmgz"];d=Max[0.25`,N[((nmax-nmin)/10)]];mins=Union[Table[FindMinimum[fn,{var,i,nmin,nmax}],{i,nmin,nmax,d}]];maxs=Union[Table[FindMinimum[-fn,{var,i,nmin,nmax}],{i,nmin,nmax,d}]];both=var/.Last[Transpose[Select[mins⋃maxs,Head[#1]==List&]]];If[both===var,both={nmin,nmax},both=N[{nmin,nmax}⋃(((Ceiling[(10^6) #1])/(10^6))&)/@both]];fnx=Function[var,fn];hghts=Sort[Transpose[{fnx/@both,both}]];fmin=Reverse[hghts[[1]]];fmax=Reverse[Last[hghts]];On[FindMinimum::"fmlim",FindMinimum::"regex",FindMinimum::"fmns",Transpose::"nmtx",FindMinimum::"fmgz"];{fmin,fmax}]

PolyMinMax[fn_,{var_,min_,max_}]:=Module[{zeros={},pts,fnx,xx},If[Exponent[fn,var]>1,zeros=Chop[var/.NSolve[(∂_var)fn==0,var]];zeros=Select[zeros,Head[#1]=!=Complex&];zeros=Select[zeros,min<#1&&#1<max&]];fnx=Function[var,fn];pts=Flatten[AppendTo[zeros,{min,max}]];pts=Sort[Transpose[Reverse[{pts,fnx/@pts}]]];pts={pts[[1]],pts[[Length[pts]]]};Reverse/@pts]

LeftApprox[fn_,{var_,min_,max_},n_:8,opts___]:=left[{0,fn},{var,min,max},n,opts]

LeftApprox[{fn1_,fn2_},{var_,min_,max_},n_:8,opts___]:=left[{fn1,fn2},{var,min,max},n,opts]

leftR[fn_,{var_,min_,max_},n_:8,color_:ColorData["Legacy", "Burlywood"]]:=leftR[{0,fn},{var,min,max},n,color]

leftR[{fn1_,fn2_},{var_,min_,max_},n_:8,color_:ColorData["Legacy", "Burlywood"]]:=Module[{xpoints,dx,i,hgts1,hgts2,llist,rlist},dx=((max-min)/n);xpoints=Table[i,{i,min,max,dx}];hgts1=Function[var,fn1]/@xpoints;hgts2=Function[var,fn2]/@xpoints;llist=Table[{Line[{{xpoints[[i]],hgts2[[i]]},{xpoints[[i+1]],hgts2[[i]]}}],Line[{{xpoints[[i+1]],hgts2[[i]]},{xpoints[[i+1]],hgts1[[i]]}}],Line[{{xpoints[[i+1]],hgts1[[i]]},{xpoints[[i]],hgts1[[i]]}}],Line[{{xpoints[[i]],hgts1[[i]]},{xpoints[[i]],hgts2[[i]]}}]},{i,n}];rlist={color,Table[Rectangle[{xpoints[[i]],hgts2[[i]]},{xpoints[[i+1]],hgts1[[i]]}],{i,1,n}]};{rlist,llist}]

left[fn_,{var_,min_,max_},n_:8,opts___]:=left[{0,fn},{var,min,max},n,opts]

left[{fn1_,fn2_},{var_,min_,max_},n_:8,opts___]:=Module[{xpoints,dx,i,hgts,ptstoo,area},ptstoo=PtsInfoToo/.{opts}/.Options[left];dx=((max-min)/n);xpoints=Drop[Table[i,{i,min,max,dx}],{n+1}];hgts=(Function[var,fn2][#1]-Function[var,fn1][#1]&)/@xpoints;area=Plus@@(hgts dx);If[ptstoo,{xpoints,area},area]]

MidpointApprox[fn_,{var_,min_,max_},n_:8,opts___]:=midpoint[{0,fn},{var,min,max},n,opts]

MidpointApprox[{fn1_,fn2_},{var_,min_,max_},n_:8,opts___]:=midpoint[{fn1,fn2},{var,min,max},n,opts]

midpointR[fn_,{var_,min_,max_},n_:8,color_:Magenta]:=midpointR[{0,fn},{var,min,max},n,color]

midpointR[{fn1_,fn2_},{var_,min_,max_},n_:8,color_:Magenta]:=Module[{xpoints,dx,i,hgts1,hgts2,llist,rlist,mdpts},dx=((max-min)/n);xpoints=Table[i,{i,min,max,dx}];mdpts=(1/2) (Drop[xpoints,{1}]+Drop[xpoints,{n+1}]);hgts1=Function[var,fn1]/@mdpts;hgts2=Function[var,fn2]/@mdpts;llist=Table[{Line[{{xpoints[[i]],hgts2[[i]]},{xpoints[[i+1]],hgts2[[i]]}}],Line[{{xpoints[[i+1]],hgts2[[i]]},{xpoints[[i+1]],hgts1[[i]]}}],Line[{{xpoints[[i+1]],hgts1[[i]]},{xpoints[[i]],hgts1[[i]]}}],Line[{{xpoints[[i]],hgts1[[i]]},{xpoints[[i]],hgts2[[i]]}}]},{i,1,n}];rlist={color,Table[Rectangle[{xpoints[[i]],hgts2[[i]]},{xpoints[[i+1]],hgts1[[i]]}],{i,1,n}]};{rlist,llist}]

midpoint[fn_,{var_,min_,max_},n_:8,opts___]:=midpoint[{0,fn},{var,min,max},n,opts]

midpoint[{fn1_,fn2_},{var_,min_,max_},n_:8,opts___]:=Module[{xpoints,dx,i,hgts,area,ptstoo},ptstoo=PtsInfoToo/.{opts}/.Options[midpoint];dx=((max-min)/n);xpoints=Table[i,{i,min,max,dx}];xpoints=(1/2) (Drop[xpoints,{1}]+Drop[xpoints,{n+1}]);hgts=(Function[var,fn2][#1]-Function[var,fn1][#1]&)/@xpoints;area=Plus@@(hgts dx);If[ptstoo,{xpoints,area},area]]

midpoint[{fn1_,fn2_},{var1_,a_,b_},{var2_,c_,d_},{n1_Integer:8,n2_Integer:8},opts___]:=Module[{xypoints,dx,dy,i,hgts,volume,ptstoo},ptstoo=PtsInfoToo/.{opts}/.Options[midpoint];dx=((b-a)/n1);dy=((d-c)/n2);xypoints=Flatten[Table[{a+(1/2) (2 i-1) dx,c+(1/2) (2 j-1) dy},{i,n1},{j,n2}],1];hgts=((fn2/.{var1->First[#1],var2->Last[#1]})-(fn1/.{var1->First[#1],var2->Last[#1]})&)/@xypoints;volume=Plus@@(hgts dx dy);If[ptstoo,{xypoints,volume},volume]]

midpoint[{fn1_,fn2_},{var1_,a_,b_},{var2_,c_,d_},n_Integer:8,opts___]:=midpoint[{fn1,fn2},{var1,a,b},{var2,c,d},{n,n},opts]

midpoint[fn_,{var1_,a_,b_},{var2_,c_,d_},n_Integer:8,opts___]:=midpoint[{0,fn},{var1,a,b},{var2,c,d},{n,n},opts]

midpoint[fn_,{var1_,a_,b_},{var2_,c_,d_},{n1_Integer:8,n2_Integer:8},opts___]:=midpoint[{0,fn},{var1,a,b},{var2,c,d},{n1,n2},opts]

nmidpoint[{fn1_,fn2_},{var1_,a_,b_},{var2_,c_,d_},{n1_Integer:8,n2_Integer:8},opts___]:=N[midpoint[{fn1,fn2},{var1,a,b},{var2,c,d},{n1,n2},opts]]

nmidpoint[{fn1_,fn2_},{var1_,a_,b_},{var2_,c_,d_},n_Integer:8,opts___]:=N[midpoint[{fn1,fn2},{var1,a,b},{var2,c,d},{n,n},opts]]

nmidpoint[fn2_,{var1_,a_,b_},{var2_,c_,d_},{n1_Integer:8,n2_Integer:8},opts___]:=N[midpoint[{0,fn2},{var1,a,b},{var2,c,d},{n1,n2},opts]]

nmidpoint[fn2_,{var1_,a_,b_},{var2_,c_,d_},n_Integer:8,opts___]:=N[midpoint[{0,fn2},{var1,a,b},{var2,c,d},{n,n},opts]]

ncircumscribed[fn_,{var_,min_,max_},n_:8,opts___]:=N[circumscribed[fn,{var,min,max},n,opts]]

ncircumscribed[{fn1_,fn2_},{var_,min_,max_},n_:8,opts___]:=N[circumscribed[{fn1,fn2},{var,min,max},n,opts]]

ninscribed[fn_,{var_,min_,max_},n_:8,opts___]:=N[inscribed[fn,{var,min,max},n,opts]]

ninscribed[{fn1_,fn2_},{var_,min_,max_},n_:8,opts___]:=N[inscribed[{fn1,fn2},{var,min,max},n,opts]]

NLeftApprox[fn_,{var_,min_,max_},n_:8,opts___]:=nleft[{0,fn},{var,min,max},n,opts]

NLeftApprox[{fn1_,fn2_},{var_,min_,max_},n_:8,opts___]:=nleft[{fn1,fn2},{var,min,max},n,opts]

nleft[fn_,{var_,min_,max_},n_:8,opts___]:=N[left[fn,{var,min,max},n,opts]]

NMidpointApprox[fn_,{var_,min_,max_},n_:8,opts___]:=nmidpoint[{0,fn},{var,min,max},n,opts]

NMidpointApprox[{fn1_,fn2_},{var_,min_,max_},n_:8,opts___]:=nmidpoint[{fn1,fn2},{var,min,max},n,opts]

nmidpoint[fn_,{var_,min_,max_},n_:8,opts___]:=N[midpoint[fn,{var,min,max},n,opts]]

nmidpoint[{fn1_,fn2_},{var_,min_,max_},n_:8,opts___]:=N[midpoint[{fn1,fn2},{var,min,max},n,opts]]

nrandom[fn_,{var_,min_,max_},n_:8,opts___]:=N[random[fn,{var,min,max},n,opts]]

nrandom[fn_,{var_,min_,max_},n_:8,opts___]:=N[random[fn,{var,min,max},n,opts]]

nrandom[{fn1_,fn2_},{var_,min_,max_},n_:8,opts___]:=N[random[{fn1,fn2},{var,min,max},n,opts]]

NRightApprox[fn_,{var_,min_,max_},n_:8,opts___]:=nright[{0,fn},{var,min,max},n,opts]

NRightApprox[{fn1_,fn2_},{var_,min_,max_},n_:8,opts___]:=nright[{fn1,fn2},{var,min,max},n,opts]

nright[fn_,{var_,min_,max_},n_:8,opts___]:=N[right[fn,{var,min,max},n,opts]]

nright[{fn1_,fn2_},{var_,min_,max_},n_:8,opts___]:=N[right[{fn1,fn2},{var,min,max},n,opts]]

NSimpsonApprox[fn_,{var_,min_,max_},n_?EvenQ,opts___]:=N[SimpsonApprox[fn,{var,min,max},n,opts]]

nsimpson[fn_,{var_,min_,max_},n_?EvenQ,opts___]:=N[simpson[fn,{var,min,max},n,opts]]

NTrapezoidApprox[fn_,{var_,min_,max_},n_:8,opts___]:=N[TrapezoidApprox[fn,{var,min,max},n,opts]]

ntrapezoid[fn_,{var_,min_,max_},n_:8,opts___]:=N[trapezoid[fn,{var,min,max},n,opts]]

Options[circumscribed]={PtsInfoToo->False};
Options[inscribed]={PtsInfoToo->False};
Options[left]={PtsInfoToo->False};
Options[midpoint]={PtsInfoToo->False};
Options[random]={PtsInfoToo->False};
Options[right]={PtsInfoToo->False};
Options[simpson]={PtsInfoToo->False};
Options[trapezoid]={PtsInfoToo->False};

randomR[fn_,{var_,min_,max_},n_:8,color_:Red]:=randomR[{0,fn},{var,min,max},n,color]

randomR[{fn1_,fn2_},{var_,min_,max_},n_:8,color_:Red]:=Module[{xpoints,dx,i,hgts1,hgts2,llist,rlist,rpts,hgts,area},dx=N[((max-min)/n)];xpoints=Table[i,{i,min,max,dx}];rpts=Table[RandomReal[{N[min+(i-1) dx],N[min+i dx]}],{i,1,n}];hgts=(Function[var,fn2][#1]-Function[var,fn1][#1]&)/@xpoints;area=N[Plus@@(hgts dx)];hgts1=Function[var,fn1]/@rpts;hgts2=Function[var,fn2]/@rpts;llist=Table[{Line[{{xpoints[[i]],hgts2[[i]]},{xpoints[[i+1]],hgts2[[i]]}}],Line[{{xpoints[[i+1]],hgts2[[i]]},{xpoints[[i+1]],hgts1[[i]]}}],Line[{{xpoints[[i+1]],hgts1[[i]]},{xpoints[[i]],hgts1[[i]]}}],Line[{{xpoints[[i]],hgts1[[i]]},{xpoints[[i]],hgts2[[i]]}}]},{i,1,n}];rlist={color,Table[Rectangle[{xpoints[[i]],hgts2[[i]]},{xpoints[[i+1]],hgts1[[i]]}],{i,1,n}]};
{area,{rlist,llist}}]

random[fn_,{var_,min_,max_},n_:8,opts___]:=random[{0,fn},{var,min,max},n,opts]

random[{fn1_,fn2_},{var_,min_,max_},n_:8,opts___]:=Module[{xpoints,dx,i,hgts,area,ptstoo},ptstoo=PtsInfoToo/.{opts}/.Options[random];dx=N[((max-min)/n)];xpoints=Table[RandomReal[{N[min+(i-1) dx],N[min+i dx]}],{i,1,n}];hgts=(Function[var,fn2][#1]-Function[var,fn1][#1]&)/@xpoints;area=Plus@@(hgts dx);If[ptstoo,{xpoints,area},area]]

RightApprox[fn_,{var_,min_,max_},n_:8,opts___]:=right[{0,fn},{var,min,max},n,opts]

RightApprox[{fn1_,fn2_},{var_,min_,max_},n_:8,opts___]:=right[{fn1,fn2},{var,min,max},n,opts]

rightR[fn_,{var_,min_,max_},n_:8,color_:Blue]:=rightR[{0,fn},{var,min,max},n,color]

rightR[{fn1_,fn2_},{var_,min_,max_},n_:8,color_:Blue]:=Module[{xpoints,dx,i,hgts1,hgts2,llist,rlist},dx=((max-min)/n);xpoints=Table[i,{i,min,max,dx}];hgts1=Function[var,fn1]/@xpoints;hgts2=Function[var,fn2]/@xpoints;llist=Table[{Line[{{xpoints[[i]],hgts2[[i+1]]},{xpoints[[i+1]],hgts2[[i+1]]}}],Line[{{xpoints[[i+1]],hgts2[[i+1]]},{xpoints[[i+1]],hgts1[[i+1]]}}],Line[{{xpoints[[i+1]],hgts1[[i+1]]},{xpoints[[i]],hgts1[[i+1]]}}],Line[{{xpoints[[i]],hgts1[[i+1]]},{xpoints[[i]],hgts2[[i+1]]}}]},{i,1,n}];rlist={color,Table[Rectangle[{xpoints[[i]],hgts2[[i+1]]},{xpoints[[i+1]],hgts1[[i+1]]}],{i,1,n}]};{rlist,llist}]

right[fn_,{var_,min_,max_},n_:8,opts___]:=right[{0,fn},{var,min,max},n,opts]

right[{fn1_,fn2_},{var_,min_,max_},n_:8,opts___]:=Module[{xpoints,dx,i,hgts,area,ptstoo},ptstoo=PtsInfoToo/.{opts}/.Options[right];dx=((max-min)/n);xpoints=Drop[Table[i,{i,min,max,dx}],{1}];hgts=(Function[var,fn2][#1]-Function[var,fn1][#1]&)/@xpoints;area=Plus@@(hgts dx);If[ptstoo,{xpoints,area},area]]

ShowApproximations2[fn_,{var_,min_,max_},type_,n_:8,opts___]:=Module[{k,rlist,llist,randlist,prt,thex,deltax,mycolors,fnx,xxx,gr1,gr2},randlist=Table[RandomReal[],{k,1,n}];deltax=((max-min)/n);fnx=Function[xxx,fn/.var->xxx];prt[0]=min;prt[m_]:=prt[m]=prt[m-1]+deltax;thex[m_]:=Switch[type,
left,prt[m-1],
right,prt[m],
midpoint,(1/2) (prt[m-1]+prt[m]),
random,prt[m-1]+randlist[[m]] deltax,
inscribed,FindAbsMinMax[fnx[var],{var,prt[m-1],prt[m]}][[1,1]],
circumscribed,FindAbsMinMax[fnx[var],{var,prt[m-1],prt[m]}][[2,1]]];mycolors=Switch[type,
left,ColorData["Legacy", "Burlywood"],
right,Blue,
midpoint,Magenta,
random,Red,
inscribed,Cyan,
circumscribed,Yellow,
trapezoid,Green,
simpson,Pink];llist=Line[Flatten[Table[{{prt[k-1],0},{prt[k-1],fnx[thex[k]]},{prt[k],fnx[thex[k]]},{prt[k],0}},{k,1,n}],1]];rlist=Table[Rectangle[{prt[k-1],0},{prt[k],fnx[thex[k]]}],{k,1,n}];gr1=Graphics[{{mycolors,rlist},{Thickness[0.005`],llist}}];gr2=Plot[fnx[var],{var,min,max},opts,PlotStyle->{{Black,Thickness[0.01`]}},PlotLabel->"("<>ToString[n]<>") Estimated area: "<>ToString[(∑_(k=1)^n)N[(prt[k]-prt[k-1]) fnx[thex[k]]]],Axes->True,DisplayFunction->Identity];Show[{gr2,gr1,gr2},opts,DisplayFunction->$DisplayFunction]]
ShowApproximations[fn_,{var_,min_,max_},type_,n_:8,opts___]:=ShowApproximations[{0,fn},{var,min,max},type,n,opts]

ShowApproximations[{fn1_,fn2_},{var1_,a_,b_},{var2_,c_,d_},type_,{n1_:8,n2_:8},opts___]:=Module[{k,rlist,llist,randlist,prt,thex,deltax,mycolors,fnx,xxx,gr1,gr2,volume,grlist,daQ,rc,finopts={opts}},daQ=DisplayArea/.{opts}/.Options[ShowApproximations];rc=Color/.{opts}/.Options[ShowApproximations];finopts=DeleteCases[finopts,Color->_];finopts=DeleteCases[finopts,DisplayArea->_];AppendTo[finopts,Axes->True];mycolors=Switch[type,
closest,Blue,
midpoint,Magenta,
random,Red,
inscribed,Cyan,
circumscribed,Yellow];If[rc=!=Black,mycolors=rc];{volume,grlist}=Switch[type,
closest,{nclosest[{fn1,fn2},{var1,a,b},{var2,c,d},{n1,n2}],closestC[fn2,{var1,a,b},{var2,c,d},{n1,n2}]},
midpoint,{nmidpoint[{fn1,fn2},{var1,a,b},{var2,c,d},{n1,n2}],midpointC[fn2,{var1,a,b},{var2,c,d},{n1,n2}]},
random,randomR[{fn1,fn2},{var,min,max},n,mycolors],
inscribed,{ninscribed[{fn1,fn2},{var,min,max},n],inscribedR[{fn1,fn2},{var,min,max},n,mycolors]},
circumscribed,{ncircumscribed[{fn1,fn2},{var,min,max},n],circumscribedR[{fn1,fn2},{var,min,max},n,mycolors]}];gr1=Show[Graphics3D[grlist],DisplayFunction->Identity,Lighting->None];gr2=Plot3D[fn2,{var1,a,b},{var2,c,d},ColorFunction->Hue,DisplayFunction->Identity];Show[{gr2,gr1,gr2},DisplayFunction->$DisplayFunction,Lighting->None,opts]]

ShowApproximations[fn_,{var1_,a_,b_},{var2_,c_,d_},type_,{n1_:8,n2_:8},opts___]:=ShowApproximations[{0,fn},{var1,a,b},{var2,c,d},type,{n1,n2},opts]

ShowSimpsonApprox[fn_,{var_,min_,max_},n_:8,opts___]:=Module[{k,rlist,llist,prt,thex,parpart,deltax,mycolors,fnx,xxx,gr1,gr2, plotrange},deltax=((max-min)/n);
plotrange=AbsoluteOptions[Plot[fn,{var, min, max}],PlotRange];
If[OddQ[n],Print["You must choose an even number of subintervals."]];fnx=Function[xxx,fn/.var->xxx];Do[prt[i]=min+i deltax,{i,0,n}];thex[m_]:=prt[m-1];llist=Table[Line[{{prt[k],0},{prt[k],fnx[prt[k]]}}],{k,0,n,2}];AppendTo[llist,Table[{Dashing[{0.01`,0.01`}],Line[{{prt[k],0},{prt[k],fnx[prt[k]]}}]},{k,1,n-1,2}]];parpart=Flatten[Table[MakeOneParabola[fn,{var,prt[i],prt[i+2]},First[plotrange]],{i,0,n-2,2}]];gr1=Graphics[{Thickness[0.005`],llist}];Show[{parpart,gr1},opts,PlotLabel->"("<>ToString[n]<>") Estimated area: "<>ToString[NSimpsonApprox[fn,{var,min,max},n]],DisplayFunction->$DisplayFunction]]
MakeOneParabola[fn_,{var_,first_,second_},opts___]:=Module[{fnx,x1,x2,x3,a,b,c,sol,gr1,gr2},fnx=Function[var,fn];x1=first;x2=((first+second)/2);x3=second;sol=NSolve[{fnx[x1]==a (x1^2)+b x1+c,fnx[x2]==a (x2^2)+b x2+c,fnx[x3]==a (x3^2)+b x3+c},{a,b,c}];{a,b,c}=Flatten[{a,b,c}/.sol];gr1=Plot[fn,{var,first,second},opts];gr2=Plot[a*(var^2)+b*var+c,{var,first,second},PlotStyle->Black,Filling-> Axis,FillingStyle->Pink,opts];{gr2,gr1}]

ShowTrapezoidApprox[fn_,{var_,min_,max_},n_:8,opts___]:=Module[{k,rlist,llist,prt,thex,deltax,mycolors,fnx,xxx,gr1,gr2},deltax=((max-min)/n);fnx=Function[var,fn];prt[0]=min;prt[m_]:=prt[m]=prt[m-1]+deltax;thex[m_]:=prt[m-1];mycolors=Green;llist=Line[Flatten[Table[{{prt[k-1],0},{prt[k-1],fnx[prt[k-1]]},{prt[k],fnx[prt[k]]},{prt[k],0}},{k,1,n}],1]];rlist=Table[Polygon[{{prt[k-1],0},{prt[k-1],fnx[prt[k-1]]},{prt[k],fnx[prt[k]]},{prt[k],0},{prt[k-1],0}}],{k,1,n}];gr1=Graphics[{{mycolors,rlist},{Thickness[0.005`],llist}}];gr2=Plot[fnx[var],{var,min,max},opts,PlotStyle->{{Black,Thickness[0.01`]}},PlotLabel->"("<>ToString[n]<>") Estimated area: "<>ToString[NTrapezoidApprox[fn,{var,min,max},n]],Axes->True,DisplayFunction->Identity];Show[{gr2,gr1,gr2},opts,DisplayFunction->$DisplayFunction]]

SimpsonApprox[fn_,{var_,min_,max_},n_?EvenQ,opts___]:=(1/3) (TrapezoidApprox[fn,{var,min,max},(n/2)]+2 MidpointApprox[fn,{var,min,max},(n/2)])

simpsonR[fn_,{var_,min_,max_},n_?EvenQ,color_:Green]:=Module[{xpoints,dx,i,hgts,llist,parpart,plotrange,rlist,gr},dx=((max-min)/n);xpoints=Table[i,{i,min,max,dx}];
plotrange=AbsoluteOptions[Plot[fn,{var, min, max}],PlotRange];hgts=Function[var,fn]/@xpoints;llist=Table[Line[{{xpoints[[i]],0},{xpoints[[i]],hgts[[i]]}}],{i,1,n+1,2}];AppendTo[llist,Table[{Dashing[{0.01`,0.01`}],Line[{{xpoints[[i]],0},{xpoints[[i]],hgts[[i]]}}]},{i,2,n,2}]];gr=Graphics[{Thickness[0.005`],llist}];parpart=Flatten[Table[MakeOneParabola[fn,{var,xpoints[[i]],xpoints[[i+2]]},First[plotrange]],{i,1,n-1,2}]];Flatten[{parpart,gr}]]

simpson[fn_,{var_,min_,max_},n_?EvenQ,opts___]:=Module[{ptstoo,xpoints,area,dx,i},ptstoo=PtsInfoToo/.{opts}/.Options[circumscribed];xpoints=If[ptstoo,dx=((max-min)/n);Table[i,{i,min,max,dx}]];area=(1/3) (trapezoid[fn,{var,min,max},(n/2)]+2 midpoint[fn,{var,min,max},(n/2)]);If[ptstoo,{xpoints,area},area]]

TrapezoidApprox[fn_,{var_,min_,max_},n_:8,opts___]:=(1/2) (LeftApprox[fn,{var,min,max},n]+RightApprox[fn,{var,min,max},n])
trapezoidR[fn_,{var_,min_,max_},n_:8,color_:Green]:=Module[{xpoints,dx,i,hgts,llist,rlist},dx=((max-min)/n);xpoints=N[Table[i,{i,min,max,dx}]];hgts=Function[var,fn]/@xpoints;llist=Table[{Line[{{xpoints[[i]],0},{xpoints[[i]],hgts[[i]]}}],Line[{{xpoints[[i]],hgts[[i]]},{xpoints[[i+1]],hgts[[i+1]]}}],Line[{{xpoints[[i+1]],hgts[[i+1]]},{xpoints[[i+1]],0}}],Line[{{xpoints[[i+1]],0},{xpoints[[i]],0}}]},{i,1,n}];rlist={color,Table[Polygon[{{xpoints[[i]],0},{xpoints[[i]],hgts[[i]]},{xpoints[[i+1]],hgts[[i+1]]},{xpoints[[i+1]],0}}],{i,1,n}]};{rlist,llist}]
trapezoid[fn_,{var_,min_,max_},n_:8,opts___]:=Module[{ptstoo,xpoints,area,dx,i},ptstoo=PtsInfoToo/.{opts}/.Options[circumscribed];xpoints=If[ptstoo,dx=((max-min)/n);Table[i,{i,min,max,dx}]];area=(1/2) (left[fn,{var,min,max},n]+right[fn,{var,min,max},n]);If[ptstoo,{xpoints,area},area]]

ShowLeftRightDifference[fn_,{var_,min_,max_},n_:8]:=Module[{increasingQ,gr1,gr2,gr3},increasingQ=((∂_var)fn/.var->((max+min)/2))>0;gr1=Graphics[leftR[fn,{var,min,max},n]];gr2=Graphics[rightR[fn,{var,min,max},n]];gr3=Plot[fn,{var,min,max},PlotLabel->"differences",PlotStyle->{{Black,Thickness[0.01`]}},Axes->True];If[increasingQ,Show[{gr2,gr1,gr3},Axes->True, AspectRatio->Automatic],Show[{gr1,gr2,gr3},Axes->True, AspectRatio->Automatic]]];

leftPts[pts_,n_, valQ_]:=Module[{rects,len,domain,a,b,dx,index,f},
len = Length[pts]-1;
If[len/n==Floor[len/n],
domain = Sort[First[Transpose[data]]];
Do[f[data[[k]]//First]=data[[k]]//Last,{k,1,Length[data]}];
a = First[domain];
b=Last[domain];
dx =(b-a)/n;
index=1;
rects = Table[Rectangle[{a+(k-1)*dx,0},{a+k*dx,f[a+(k-1)*dx]}],{k,index,n}];
lines = Table[Line[{{a+(k-1)*dx,0},{a+(k-1)*dx,f[a+(k-1)*dx]},{a+k*dx,f[a+(k-1)*dx]},{a+k*dx,0}}],{k,index,n}];
area = Sum[dx *f[a+(k-1)*dx],{k,index,n}];
rects = {Thick,ColorData["Legacy", "Burlywood"],rects,Black,lines};
If[valQ,ListPlot[data,PlotStyle->{Blue,PointSize[Large]},AxesOrigin->{0,0},Prolog->rects,PlotLabel->ToString[area]],
ListPlot[data,PlotStyle->{Blue,PointSize[Large]},AxesOrigin->{0,0},Prolog->rects]]
]
]
rightPts[pts_,n_, valQ_]:=Module[{rects,len,domain,a,b,dx,index,f},
len = Length[pts]-1;
If[len/n==Floor[len/n],
domain = Sort[First[Transpose[data]]];
Do[f[data[[k]]//First]=data[[k]]//Last,{k,1,Length[data]}];
a = First[domain];
b=Last[domain];
dx =(b-a)/n;
index=1;
rects = Table[Rectangle[{a+(k-1)*dx,0},{a+k*dx,f[a+k*dx]}],{k,index,n}];
lines = Table[Line[{{a+(k-1)*dx,0},{a+(k-1)*dx,f[a+k*dx]},{a+k*dx,f[a+k*dx]},{a+k*dx,0}}],{k,index,n}];
area = Sum[dx *f[a+k*dx],{k,index,n}];
rects = {Thick,Cyan,rects,Black,lines};
If[valQ,ListPlot[data,PlotStyle->{Blue,PointSize[Large]},AxesOrigin->{0,0},Prolog->rects,PlotLabel->ToString[area]],
ListPlot[data,PlotStyle->{Blue,PointSize[Large]},AxesOrigin->{0,0},Prolog->rects]]
]
]
midPts[pts_,n_, valQ_]:=Module[{rects,len,domain,a,b,dx,index,f},
len = Length[pts]-1;
If[len/n==Floor[len/n]&&EvenQ[len/n],
domain = Sort[First[Transpose[data]]];
Do[f[data[[k]]//First]=data[[k]]//Last,{k,1,Length[data]}];
a = First[domain];
b=Last[domain];
dx =(b-a)/n;
index=1;
rects = Table[Rectangle[{a+(k-1)*dx,0},{a+k*dx,f[a+(2k-1)/2*dx]}],{k,index,n}];
lines = Table[Line[{{a+(k-1)*dx,0},{a+(k-1)*dx,f[a+(2k-1)/2*dx]},{a+k*dx,f[a+(2k-1)/2*dx]},{a+k*dx,0}}],{k,index,n}];
area = Sum[dx *f[a+(2k-1)/2*dx],{k,index,n}];
rects = {Thick,Magenta,rects,Black,lines};
If[valQ,ListPlot[data,PlotStyle->{Blue,PointSize[Large]},AxesOrigin->{0,0},Prolog->rects,PlotLabel->ToString[area]],
ListPlot[data,PlotStyle->{Blue,PointSize[Large]},AxesOrigin->{0,0},Prolog->rects]]
]
]

Options[midpoint]={PtsInfoToo->False};
midpoint[{fn1_,fn2_},{var1_,a_,b_},{var2_,c_,d_},{n1_Integer:8,n2_Integer:8},opts___]:=Module[{xypoints,dx,dy,i,hgts,volume,ptstoo},ptstoo=PtsInfoToo/.{opts}/.Options[midpoint];dx=((b-a)/n1);dy=((d-c)/n2);xypoints=Flatten[Table[{a+(1/2) (2 i-1) dx,c+(1/2) (2 j-1) dy},{i,n1},{j,n2}],1];hgts=((fn2/.{var1->First[#1],var2->Last[#1]})-(fn1/.{var1->First[#1],var2->Last[#1]})&)/@xypoints;volume=Plus@@(hgts dx dy);If[ptstoo,{xypoints,volume},volume]]
midpoint[{fn1_,fn2_},{var1_,a_,b_},{var2_,c_,d_},n_Integer:8,opts___]:=midpoint[{fn1,fn2},{var1,a,b},{var2,c,d},{n,n},opts]
midpoint[fn_,{var1_,a_,b_},{var2_,c_,d_},n_Integer:8,opts___]:=midpoint[{0,fn},{var1,a,b},{var2,c,d},{n,n},opts]
midpoint[fn_,{var1_,a_,b_},{var2_,c_,d_},{n1_Integer:8,n2_Integer:8},opts___]:=midpoint[{0,fn},{var1,a,b},{var2,c,d},{n1,n2},opts]
nmidpoint[{fn1_,fn2_},{var1_,a_,b_},{var2_,c_,d_},{n1_Integer:8,n2_Integer:8},opts___]:=N[midpoint[{fn1,fn2},{var1,a,b},{var2,c,d},{n1,n2},opts]]
nmidpoint[{fn1_,fn2_},{var1_,a_,b_},{var2_,c_,d_},n_Integer:8,opts___]:=N[midpoint[{fn1,fn2},{var1,a,b},{var2,c,d},{n,n},opts]]
nmidpoint[fn2_,{var1_,a_,b_},{var2_,c_,d_},{n1_Integer:8,n2_Integer:8},opts___]:=N[midpoint[{0,fn2},{var1,a,b},{var2,c,d},{n1,n2},opts]]
nmidpoint[fn2_,{var1_,a_,b_},{var2_,c_,d_},n_Integer:8,opts___]:=N[midpoint[{0,fn2},{var1,a,b},{var2,c,d},{n,n},opts]]
Options[closest]={PtsInfoToo->False};
closest[{fn1_,fn2_},{var1_,a_,b_},{var2_,c_,d_},{n1_Integer:8,n2_Integer:8},opts___]:=Module[{xypoints,dx,dy,i,hgts,volume,ptstoo},ptstoo=PtsInfoToo/.{opts}/.Options[closest];dx=((b-a)/n1);dy=((d-c)/n2);xypoints=Flatten[Table[{a+(i-1) dx,c+(j-1) dy},{i,n1},{j,n2}],1];hgts=((fn2/.{var1->First[#1],var2->Last[#1]})-(fn1/.{var1->First[#1],var2->Last[#1]})&)/@xypoints;volume=Plus@@(hgts dx dy);If[ptstoo,{xypoints,volume},volume]]
closest[{fn1_,fn2_},{var1_,a_,b_},{var2_,c_,d_},n_Integer:8,opts___]:=closest[{fn1,fn2},{var1,a,b},{var2,c,d},{n,n},opts]
closest[fn_,{var1_,a_,b_},{var2_,c_,d_},n_Integer:8,opts___]:=closest[{0,fn},{var1,a,b},{var2,c,d},{n,n},opts]
closest[fn_,{var1_,a_,b_},{var2_,c_,d_},{n1_Integer:8,n2_Integer:8},opts___]:=closest[{0,fn},{var1,a,b},{var2,c,d},{n1,n2},opts]
nclosest[{fn1_,fn2_},{var1_,a_,b_},{var2_,c_,d_},{n1_Integer:8,n2_Integer:8},opts___]:=N[closest[{fn1,fn2},{var1,a,b},{var2,c,d},{n1,n2},opts]]
nclosest[{fn1_,fn2_},{var1_,a_,b_},{var2_,c_,d_},n_Integer:8,opts___]:=N[closest[{fn1,fn2},{var1,a,b},{var2,c,d},{n,n},opts]]
nclosest[fn2_,{var1_,a_,b_},{var2_,c_,d_},{n1_Integer:8,n2_Integer:8},opts___]:=N[closest[{0,fn2},{var1,a,b},{var2,c,d},{n1,n2},opts]]
nclosest[fn2_,{var1_,a_,b_},{var2_,c_,d_},n_Integer:8,opts___]:=N[closest[{0,fn2},{var1,a,b},{var2,c,d},{n,n},opts]]
midpointC[f_,{var1_,a_,b_},{var2_,c_,d_},n_Integer,opts___]:=midpointC[f,{var1,a,b},{var2,c,d},{n,n},opts]
midpointC[f_,{var1_,a_,b_},{var2_,c_,d_},{n1_Integer,n2_Integer},opts___]:=Module[{dx,dy,lowers,uppers,cuboids},dx=((b-a)/n1);dy=((d-c)/n2);lowers=Flatten[Table[{a+(i-1) dx,c+(j-1) dy,0},{i,n1},{j,n2}],1];uppers=Flatten[Table[{a+i dx,c+j dy,f/.{var1->a+(1/2) (2 i-1) dx,var2->c+(1/2) (2 j-1) dy}},{i,n1},{j,n2}],1];cuboids=Transpose[{lowers,uppers}];Prepend[(Cuboid[First[#1],Last[#1]]&)/@cuboids,RGBColor[1,1,1]]]
closestC[f_,{var1_,a_,b_},{var2_,c_,d_},n_Integer,opts___]:=closestC[f,{var1,a,b},{var2,c,d},{n,n},opts]
closestC[f_,{var1_,a_,b_},{var2_,c_,d_},{n1_Integer,n2_Integer},opts___]:=Module[{dx,dy,lowers,uppers,cuboids},dx=((b-a)/n1);dy=((d-c)/n2);lowers=Flatten[Table[{a+(i-1) dx,c+(j-1) dy,0},{i,n1},{j,n2}],1];uppers=Flatten[Table[{a+i dx,c+j dy,f/.{var1->a+(i-1) dx,var2->c+(j-1) dy}},{i,n1},{j,n2}],1];cuboids=Transpose[{lowers,uppers}];Prepend[(Cuboid[First[#1],Last[#1]]&)/@cuboids,RGBColor[1,1,1]]]
ShowApproximations[{fn1_,fn2_},{var1_,a_,b_},{var2_,c_,d_},type_,{n1_:8,n2_:8},opts___]:=Module[{k,rlist,llist,randlist,prt,thex,deltax,mycolors,fnx,xxx,gr1,gr2,volume,grlist,daQ,rc,finopts={opts}},daQ=DisplayArea/.{opts}/.Options[ShowApproximations];rc=Color/.{opts}/.Options[ShowApproximations];finopts=DeleteCases[finopts,Color->_];finopts=DeleteCases[finopts,DisplayArea->_];AppendTo[finopts,Axes->True];mycolors=Switch[type,
closest,Blue,
midpoint,Magenta,
random,Red,
inscribed,Cyan,
circumscribed,Yellow];If[rc=!=Black,mycolors=rc];{volume,grlist}=Switch[type,
closest,{nclosest[{fn1,fn2},{var1,a,b},{var2,c,d},{n1,n2}],closestC[fn2,{var1,a,b},{var2,c,d},{n1,n2}]},
midpoint,{nmidpoint[{fn1,fn2},{var1,a,b},{var2,c,d},{n1,n2}],midpointC[fn2,{var1,a,b},{var2,c,d},{n1,n2}]},
random,randomR[{fn1,fn2},{var,min,max},n,mycolors],
inscribed,{ninscribed[{fn1,fn2},{var,min,max},n],inscribedR[{fn1,fn2},{var,min,max},n,mycolors]},
circumscribed,{ncircumscribed[{fn1,fn2},{var,min,max},n],circumscribedR[{fn1,fn2},{var,min,max},n,mycolors]}];gr1=Show[Graphics3D[grlist],DisplayFunction->Identity,Lighting->None];gr2=Plot3D[fn2,{var1,a,b},{var2,c,d},ColorFunction->Hue,DisplayFunction->Identity];Show[{gr2,gr1,gr2},DisplayFunction->$DisplayFunction,Lighting->None,opts]]

ShowApproximations[fn_,{var1_,a_,b_},{var2_,c_,d_},type_,{n1_:8,n2_:8},opts___]:=ShowApproximations[{0,fn},{var1,a,b},{var2,c,d},type,{n1,n2},opts]
(*{1->{{2}}*)

ShowMinMax[fn_,{var_,min_,max_},opts___]:=Module[{gr,fmin,fmax,pts,fnx,hghts,xx},gr=Plot[fn,{var,min,max},Filling->Axis,FillingStyle-> Yellow,opts,DisplayFunction->Identity];fnx=Function[var,fn];pts=(xx/.NSolve[(∂_xx)fnx[xx]==0,xx])⋃{min,max};hghts=Sort[Transpose[{fnx/@pts,pts}]];fmin=hghts[[1,1]];fmax=Last[hghts][[1]];Show[{gr,Graphics[{{Blue,Line[{{min,fmin},{max,fmin}}]},{Black,Line[{{min,fmax},{max,fmax}}]}}]},DisplayFunction->$DisplayFunction,Axes->True]]

ShowGuessedArea[fn_,{var_,min_,max_},h_,opts___]:=Module[{gr,act,ar},ar=(max-min) h;act=NIntegrate[fn,{var,min,max}];mess=If[act<ar,StringJoin["Needs to be lower."],StringJoin["Needs to be higher."]];Plot[{h,fn},{var,min,max},PlotLabel->Style["height = "<>ToString[N[h]]<>"\nguessed area = "<>ToString[N[ar]]<>" "<>mess,Medium],Filling->{1->{Axis,Green},1->{{2},Yellow}},PlotStyle->{{Blue,Thickness[0.01`]},{Black,Thickness[0.01`]}},opts]]

Derivative-related

ShowSecant[fn_,{var_,min_,max_},fixed_,second_,opts___]:=Module[{pts,y,pt1,pt2},
pt1={fixed,fn/.var->fixed};
pt2={second,fn/.var->second};
pts={Blue,PointSize[0.025],Point[pt1],Magenta,Point[pt2]};
y=pt1[[2]]+N[slope[pt1,pt2]]*(var-fixed);
Plot[{fn,y},{var,min,max},opts,PlotStyle->{Black,Magenta},Epilog->pts]]

ShowSecantTan[fn_,{var_,min_,max_},fixed_List,second_List,opts___]:=Module[{pts,y,tan,fp,slp},pts={{Blue,PointSize[0.03],Point[fixed]},{Magenta,PointSize[0.025],Point[second]}};
slp=slope[fixed,second]//N;
fp=(D[fn,var]/.var->fixed[[1]])//N;
y=fixed[[2]]+slp*(var-fixed[[1]]);
tan=fixed[[2]]+fp*(var-fixed[[1]]);
Plot[{fn,y,tan},{var,min,max},opts,PlotStyle->{Black,Magenta,Blue},Epilog->pts,PlotLabel->StringJoin[" secant slope = ",ToString[slp],"\ntangent slope = ?"(*,ToString[fp]*)]]]

ShowSecantTan[fn_,{var_,min_,max_},fixed_,second_,opts___]:=
ShowSecantTan[fn,{var,min,max},{fixed,fn/.var->fixed},{second,fn/.var->second},opts]

slope[pt1_,pt2_]:=(pt2[[2]]-pt1[[2]])/(pt2[[1]]-pt1[[1]])

DifferenceQuotient[fn_,{var_,point_:a},delta_:h]:=(((fn/.var-> (point+delta)) -(fn/.var-> point))/delta)

Options[ShowSecantToTan]={Side->Either};

ShowSecantToTan[fn_,{var_,min_,max_},xpt_,n_Integer:7,opts___]:=Module[{pts,y,tan,d,fnx,pt,len,ok,side,mopts, lpts, rpts},
ok=(min≤xpt)&&(xpt≤max);
side=Side/.{opts}/.Options[ShowSecantToTan];
If[ok,lpts =closePoints[fn,var,xpt,min,n];
rpts =closePoints[fn,var,xpt,max,n];
If[Abs[min-xpt]>Abs[max-xpt],pts=lpts,pts=rpts];
If[side===Left,pts=lpts];
If[side===Right,pts=rpts];
len=Length[pts];
pt={xpt,(fn/.var->xpt)};
(*{mopts}={opts}/.{a___,Side->b_,c___}->{a,c};*)
mopts=Select[{opts},(#[[1]]=!=Side)&];
mopts=ReplacePart[mopts,Sequence,{0}];
Manipulate[ShowSecantTan[fn,{var,min,max},pt,pts[[i]],mopts],(*{i,Table[k,{k,1,len}*){i,1,len,1}]]]

(* fix the above for multiple options, both graphics and Side *)

closePoints[expr_,var_,focpt_,strtpt_,n_,fn_:none]:=Module[{pts,d=strtpt-focpt,grow,vals,y},If[fn===none,grow[y_]:=(1/2)^y,grow=fn];
pts=Table[focpt+grow[i]*d,{i,0,n-1}];
vals=Map[Function[var,expr],pts];
Transpose[{pts,vals}]]

(* plotwithtangentswork[fn_,{var_,min_,max_},pts_List,opts___]:=Module[{gr,lines={},n,len,pt,slp,dots,slopes={},newopts},n=Length[pts];
len=N[(max-min)/5];
dots=Table[Point[{pts[[i]],(fn/.var->pts[[i]])}],{i,1,n}];
lines=Table[pt=pts[[i]];
slp=D[fn,var]/.var->pt;
AppendTo[slopes,slp];
If[i<n,len=N[pts[[i+1]]-pt]/2];
Line[{{pt-3/2len,slp(-3/2len)+(fn/.var->pt)},{pt+3/2len,slp(3/2len)+(fn/.var->pt)}}],{i,1,n}];
gr=Plot[fn,{var,min,max},opts,PlotStyle->Red,DisplayFunction->Identity,Epilog->{{PointSize[0.03],Magenta,dots},{Blue,lines}}];
{gr,Transpose[{pts,slopes}]//N}] *)

plotwithtangentswork[fn_,{var_,min_,max_},pts_List,opts___]:=Module[{gr,lines={},n,len,pt,slp,dots,dots2,slopes,fpts,newopts,ar},n=Length[pts];
(* len=N[(max-min)/n]; *)
dots=Map[Point,fpts=Transpose[{pts,fn/.var->pts}]];
slopes =Map[(D[fn,var]/.var->#)&,pts];
gr=Plot[fn,{var,min,max},opts,PlotStyle->{Thickness[0.012],Red},DisplayFunction->Identity];
ar=AspectRatio/.FullOptions[gr];
lines=makeLines[{min,max},fpts,slopes];
(* Table[Line[midEndpoints[fpts[[i]],slopes[[i]],len(slopes[[i]]+1)/ar,ar]],{i,1,n}];*)
gr=Plot[fn,{var,min,max},opts,PlotStyle->{Thickness[0.012],Red},DisplayFunction->Identity,Epilog->{{PointSize[0.03],Magenta,dots},{Blue,Thickness[0.008],lines}}];
{gr,Transpose[{pts,slopes}]//N}]

plotwithtangentswork[fn_,{var_,min_,max_},n_Integer?Positive,opts___]:=plotwithtangentswork[fn,{var,min,max},Table[min+i (max-min)/n-(max-min)/(2n),{i,1,n}],opts]

(* midEndpoints[{x0_,y0_},m_,len_,ar_]:=Module[{x1,y1,rad= N[Sqrt[(len^2)/(1+(m^2))]],ends},
x1 =(x0-rad);
y1 :=( m(x1-x0)+y0);
ends = {{x1,y1}};
x1= (x0+rad);
Append[ends,{x1,y1}]
]
 dropped  - couldn't get aspect ratio to be reflected *)

makeLines[{min_,max_},fpts_List,slopes_List]:=Module[{nxpts},
nxpts = Join[{min},Table[((fpts[[k,1]]+fpts[[k+1,1]])/2),{k,1,Length[fpts]-1}],{max}];
Table[Line[{{nxpts[[k]],slopes[[k]](nxpts[[k]]-fpts[[k,1]])+fpts[[k,2]]},{nxpts[[k+1]],slopes[[k]](nxpts[[k+1]]-fpts[[k,1]])+fpts[[k,2]]}}],{k,Length[fpts]}]
]

Options[PlotWithTangents]={ShowSlopes->True};

PlotWithTangents[fn_,{var_,min_,max_},pts_List,opts___]:=Module[{returned,showslopes=True,gr,finopts={opts}},showslopes=ShowSlopes/.{opts}/.Options[PlotWithTangents];
finopts=Select[{opts},(#[[1]]=!=ShowSlopes)&];
finopts=ReplacePart[finopts,Sequence,{0}];
returned=plotwithtangentswork[fn,{var,min,max},pts,finopts];
If[showslopes,Print[returned[[2]]]];
gr=Show[returned[[1]],DisplayFunction->$DisplayFunction]
]

PlotWithTangents[fn_,{var_,min_,max_},n_Integer?Positive,opts___]:=PlotWithTangents[fn,{var,min,max},Table[min+i (max-min)/n-(max-min)/(2n),{i,1,n}],opts]

parametricplotwithtangentswork[{fx_,fy_},{var_,min_,max_},pts_List,opts___]:=Module[{gr,n,dots,slopes,lines,fpts},
n=Length[pts];
fpts=Table[{fx/.var->pts[[k]],fy/.var->pts[[k]]},{k,n}];
dots=Map[Point,fpts];
slopes = Table[((D[fy,var]/.var->pts[[k]])/(D[fx,var]/.var->pts[[k]])),{k,n}];
lines=makeLines[{min,max},fpts,slopes];
gr = ParametricPlot[{fx,fy},{var, min,max},opts,PlotStyle->{Thickness[0.012],Red},DisplayFunction->Identity,Epilog->{{PointSize[0.03],Magenta,dots},{Blue,Thickness[0.008],lines}}];
{gr,Transpose[{pts,slopes}]//N}
]

parametricplotwithtangentswork[{fx_,fy_},{var_,min_,max_},n_Integer?Positive,opts___]:=parametricplotwithtangentswork[{fx,fy},{var,min,max},Table[min+i (max-min)/n-(max-min)/(2n),{i,n}],opts]

Options[ParametricPlotWithTangents]={ShowSlopes->True};

ParametricPlotWithTangents[{fx_,fy_},{var_,min_,max_},pts_List,opts___]:=Module[{returned,showslopes=True,gr,finopts={opts}},showslopes=ShowSlopes/.{opts}/.Options[ParametricPlotWithTangents];
finopts=Select[{opts},(#[[1]]=!=ShowSlopes)&];
finopts=ReplacePart[finopts,Sequence,{0}];
returned=parametricplotwithtangentswork[{fx,fy},{var,min,max},pts,finopts];
gr=Show[returned[[1]],DisplayFunction->$DisplayFunction];
If[showslopes,gr;returned[[2]],gr]]

ParametricPlotWithTangents[{fx_,fy_},{var_,min_,max_},n_Integer?Positive,opts___]:=ParametricPlotWithTangents[{fx,fy},{var,min,max},Table[min+i (max-min)/n-(max-min)/(2n),{i,1,n}],opts]

implicitplotwithtangentswork[eq_,{var_,min_,max_},{vary_,miny_,maxy_},pts_List,opts___]:=Module[{gr,lines,vertpts,posVert,len,npts,slp,dots,slopes,nslopes,hgt,newopts,othervar,pr},Off[Power::infy];
dots = Map[Point,pts];
{othervar}=Complement[Cases[eq,_Symbol,Infinity],{var}];
slopes = (Dt[othervar,var]/.Solve[Dt[eq,x],Dt[othervar,var]])/.Map[{var->First[#],othervar->Last[#]}&,pts]//Flatten;
On[Power::infy];
len = (max-min)/Length[pts//Transpose//First//Union];gr = ContourPlot[eq,{var,min,max},{vary,miny,maxy},opts,PlotStyle->{Thickness[0.012],Red},DisplayFunction->Identity];
pr = PlotRange/.FullOptions[gr];
hgt = pr[[2,2]]-pr[[2,1]];
posVert = Position[slopes,ComplexInfinity];
vertpts = pts[[posVert//Flatten]];
npts = Delete[pts,posVert];
nslopes = Delete[slopes,posVert];
lines = Table[Line[{{vertpts[[k,1]],vertpts[[k,2]]+.25 hgt},{vertpts[[k,1]],vertpts[[k,2]]-.25 hgt}}],{k,Length[posVert]}];
lines=Join[lines,Table[Line[{{npts[[k,1]]-len,npts[[k,2]]+nslopes[[k]](-len)},{npts[[k,1]]+len,npts[[k,2]]+nslopes[[k]](len)}}],{k,Length[npts]}]];
gr = ContourPlot[eq,{var,min,max},{vary,miny,maxy},opts,PlotStyle->{Thickness[0.012],Red},Epilog->{{PointSize[0.03],Magenta,dots},{Blue,Thickness[0.008],lines}}];
List[gr,Transpose[{pts,slopes}]//N]]

Options[ContourPlotWithTangents]={ShowSlopes->True};

ContourPlotWithTangents[eq_,{var_,min_,max_},{vary_,miny_,maxy_},pts_List?MatrixQ,opts___]:=Module[{returned,showslopes=True,gr,finopts={opts}},showslopes=ShowSlopes/.{opts}/.Options[ContourPlotWithTangents];
finopts=Select[{opts},(#[[1]]=!=ShowSlopes)&];
finopts=ReplacePart[finopts,Sequence,{0}];
returned=implicitplotwithtangentswork[eq,{var,min,max},{vary,miny,maxy},pts,finopts];
gr=Show[returned[[1]],DisplayFunction->$DisplayFunction];
gr=returned[[1]];
If[showslopes,Print[returned[[2]]];gr,gr]]

ShowFunctionAndDerivative[fn_,{var_,min_,max_},pts_,opts___]:=Module[{gr1,gr2,theout},{gr1,theout}=plotwithtangentswork[fn,{var,min,max},pts,opts];gr2=ListPlot[theout,Joined->True,PlotStyle->Blue,Epilog->{Magenta,PointSize[0.03`],Point/@theout},DisplayFunction->Identity,PlotRange->{{min,max},Automatic}];
Print[theout];
Show[GraphicsGrid[{{gr1},{gr2}}]]]

ShowFunctionAndDerivative[fn_,{var_,min_,max_},n_Integer?Positive,opts___]:=ShowFunctionAndDerivative[fn,{var,min,max},Table[min+((i (max-min))/n)-((max-min)/(2 n)),{i,1,n}],opts]
ShowFunctionAndDerivativeTogether[fn_,{var_,min_,max_},pts_,opts___]:=Module[{gr1,gr2,theout},{gr1,theout}=plotwithtangentswork[fn,{var,min,max},pts,opts];gr2=ListPlot[theout,Joined->True,PlotStyle->Blue,Epilog->{Magenta,PointSize[0.03`],Point/@theout},DisplayFunction->Identity];
Print[theout];
Show[{gr2,gr1},DisplayFunction->$DisplayFunction]]

ShowFunctionAndDerivativeTogether[fn_,{var_,min_,max_},n_Integer?Positive,opts___]:=ShowFunctionAndDerivativeTogether[fn,{var,min,max},Table[min+((i (max-min))/n)-((max-min)/(2 n)),{i,1,n}],opts]
MarkIntervals[fn_,{var_,min_,max_},opts___]:=Module[{ipts,spts,ilns,slns,f,x},Clear[x];f=Function[var,fn];ipts=var/.NSolve[f''[var]==0,var];ilns=Table[Line[{{ipts[[i]],-250},{ipts[[i]],250}}],{i,Length[ipts]}];spts=var/.NSolve[f'[var]==0,var];slns=Table[Line[{{spts[[i]],-250},{spts[[i]],250}}],{i,Length[spts]}];Plot[{f[var],f'[var],f''[var]},{var,min,max},Epilog->{{Magenta,ilns},{Cyan,slns}}]]
Tangraph[fn_,{var_,min_,max_},c_,opts___]:=Module[{f,x,gr1,gr2},f=Function[var,fn];gr1=Plot[{f[var],f[c]+f'[c] (var-c)},{var,min,max},opts,PlotStyle->{Red,Blue},Epilog->{{Magenta,PointSize[0.03`],Point[{c,f[c]}]},{Blue,PointSize[0.03`],Point[{c,f'[c]}]}},PlotLabel->"at "<>ToString[N[c]]<>", slope = "<>ToString[N[f'[c]]],DisplayFunction->Identity];gr2=Plot[f'[var],{var,min,c},opts,PlotStyle->Magenta,DisplayFunction->Identity];Show[{gr1,gr2},DisplayFunction->$DisplayFunction]]

Limit-related

LimitTableIllustration[exp_,{var_,target_,initdx_:0.1`},steps_:3]:=Module[{inputs,k,outputs,temp},inputs=Sort[Flatten[Table[{target-initdx (10^(-k)),target+initdx (10^(-k))},{k,0,steps-1}]]];outputs=exp/.var->inputs;PrependTo[inputs,var];inputs=Insert[inputs,target,-steps-1];PrependTo[outputs,exp];outputs=Insert[outputs,"   ",-steps-1];temp=Map[ToBoxes,Table[{inputs,outputs}],{2}];temp=temp/."\"   \""->ButtonBox["   ",Background->RGBColor[1,1,0]];DisplayForm[FrameBox[GridBox[temp,RowLines->True,ColumnLines->True]]]]
ShowEpsilonDelta[fn_,{var_,min_,max_},a_,eps_,opts___]:=Module[{L,epLines,delLines,lowera,uppera,delta,crosslines,epline,deline,slp,slpfac,h,d="δ",e="ϵ"},L=Limit[fn,var->a];slp=(∂_var)fn/.var->a;slpfac=If[slp>0,1,-1];crosslines={Blue,Line[{{min,L},{max,L}}],Line[{{a,-250},{a,250}}]};epLines={Green,Line[{{min,L-eps},{max,L-eps}}],Line[{{min,L+eps},{max,L+eps}}]};lowera=var/.FindRoot[fn==L-eps,{var,a}];uppera=var/.FindRoot[fn==L+eps,{var,a}];delta=Min[Abs[a-lowera],Abs[a-uppera]];delLines={Cyan,Line[{{a-delta,-250},{a-delta,250}}],Line[{{a+delta,-250},{a+delta,250}}]};deline={Black,Line[{{a,h=fn/.var->a+slpfac 4 delta},{a+delta,h}}],Text[d,{a+1.2` delta,h},{-1,0}],Line[{{a,h=fn/.var->a-slpfac 4 delta},{a-delta,h}}],Text[d,{a-1.2` delta,h},{1,0}]};epline={Black,Line[{{a+2 delta,L},{a+2 delta,L+eps}}],Text[e,{a+2.2` delta,L+(eps/2)},{-1,0}],Line[{{a-2 delta,L},{a-2 delta,L-eps}}],Text[e,{a-2.2` delta,L-(eps/2)},{1,0}]};Plot[fn,{var,min,max},opts,Epilog->{crosslines,epLines,epline,deline,delLines},opts]]

Misc

Options[RandomPolynomial]={Indeterminate->x,RangeMin->-8,RangeMax->8,CoefficientRing->Integer};

PlotTaylorPolynomial[fn_,{var_,varmin_,varmax_},{base_, n_},options___]:=Module[{func,m=2},
func=Normal[Series[fn,{var,base,n}]];
Plot[{fn, func,},{var,varmin,varmax},PlotStyle->{{AbsoluteThickness[3],RGBColor[0,0,0]},{AbsoluteThickness[3],RGBColor[0,0,1]}},options]]

PlotTaylorPolynomials[fn_,{var_,varmin_,varmax_},{base_, n_},options___]:=Module[{funcs},
funcs=Table[Normal[Series[fn,{var,base,k}]],{k,0,n}];
Print[funcs];
Plot[Evaluate[Append[funcs,fn]],{var,varmin,varmax},PlotStyle->Append[Table[{AbsoluteThickness[3],Hue[(k/(n+1))]},{k,n+1}],{AbsoluteThickness[4],RGBColor[0,0,0]}],options]]

TaylorPolynomial[fn_,{var_,varmin_,varmax_},{base_, n_}]:=Normal[Series[fn,{var,base,n}]]

ReturnRandom[Integer,min_,max_]:=RandomInteger[{min,max}]
ReturnRandom[Real,min_,max_]:=RandomReal[{min,max}]
ReturnRandom[Z[k_],min_,max_]:=RandomInteger[{0,k-1}]

RandomPolynomial[degree_Integer?Positive,opts___?OptionQ]:=Module[{coeff,min,max, type,ind,j},
min=RangeMin/.Flatten[{opts,Options[RandomPolynomial]}];
max=RangeMax/.Flatten[{opts,Options[RandomPolynomial]}];
type=CoefficientRing/.Flatten[{opts,Options[RandomPolynomial]}];
ind=Indeterminate/.Flatten[{opts,Options[RandomPolynomial]}];
coeff = Table[ReturnRandom[type,min,max],{degree+1}];
If[First[coeff]==0,coeff[[1]]=1];
Reverse[coeff].Table[(ind^j),{j,0,degree}]]

RainbowDisplay[{min_,max_}]:=Show[Graphics[Table[{Hue[(j/(max-min+1))],Rectangle[{min-1+j,0},{min+j,1}]},{j,min,max}]]]
SomeFigures[]:=Show[Graphics[{{Yellow,Polygon[{{0,0},{0,3},{2,0}}]},{Text["1",{0.75`,1}]},{Green,Rectangle[{1,2},{3,4}]},{Text["2",{2,3}]},{Cyan,Polygon[{{-1,5},{2.5`,5},{2,7},{0,7}}]},{Text["3",{0.75`,6}]},{Red,Polygon[{{2.25`,4.25`},{2.25`,4.75`},{3.25`,4.75`},{3.25`,6},{3.75`,6},{3.75`,4.75`},{4.75`,4.75`},{4.75`,4.25`},{3.75`,4.25`},{3.75`,0.5`},{3.25`,0.5`},{3.25`,4.25`}}]},{Text["4",{3.5`,4.5`}]},{Magenta,Rectangle[{4.5`,0},{7.5`,1}]},{Text["5",{6,0.5`}]},{Blue,Disk[{5.5`,2},2,{0,(π/3)}]},{Text["6",{6.8`,2.5`}]},{Black,Line[{{8,1.5`},{8,6.5`},{9,5.5`},{9,0.5`},{8,1.5`}}]},{Text["7",{8.5`,3.5`}]},{Orange,Polygon[{{7.`,6},{6.5`,6.866`},{5.5`,6.866`},{5.`,6},{5.5`,5.133`},{6.5`,5.133`}}]},{Text["8",{6,6}]}}],AspectRatio->Automatic]

FindZeros[f_,var_]:= Cases[var /.NSolve[f==0,var],_Real]

PlotIncreasingDecreasing[f_,{var_,min_,max_},opts___]:=Module[{fpzeros = FindZeros[D[f,var],var],eps = 0.01 ,testvals,intervals,graphs,inc},
testvals = Join[{First[fpzeros]-eps},eps+fpzeros];
inc =Positive[D[f,var]/.{var ->testvals}];
intervals = Join[{{min,First[fpzeros]}},Table[{fpzeros[[k]],fpzeros[[k+1]]},{k,1,Length[fpzeros]-1}],{{Last[fpzeros],max}}];
Do[graphs[k] = Plot[f,{var,intervals[[k,1]],intervals[[k,2]]},opts,PlotStyle->{Thickness[0.01],If[inc[[k]],RGBColor[1,0,1],RGBColor[.4,.4,.4]]},DisplayFunction->Identity],{k,1,Length[fpzeros]+1}];
Show[Table[graphs[k],{k,Length[fpzeros]+1}],DisplayFunction->$DisplayFunction]
]

PlotIncreasingDecreasing[{f_,g_},{var_,min_,max_},opts___]:=Module[{fpzeros = FindZeros[D[f,var],var],eps = 0.01 ,testvals,intervals,graphs,inc, ggr},
testvals = Join[{First[fpzeros]-eps},eps+fpzeros];
inc =Positive[D[f,var]/.{var ->testvals}];
intervals = Join[{{min,First[fpzeros]}},Table[{fpzeros[[k]],fpzeros[[k+1]]},{k,1,Length[fpzeros]-1}],{{Last[fpzeros],max}}];
ggr=Plot[g,{var, min,max},PlotStyle-> {Thickness[0.01],RGBColor[0,1,0]},DisplayFunction->Identity];
Do[graphs[k] = Plot[f,{var,intervals[[k,1]],intervals[[k,2]]},opts,PlotStyle->{Thickness[0.01],If[inc[[k]],RGBColor[1,0,1],RGBColor[.4,.4,.4]]},DisplayFunction->Identity],{k,1,Length[fpzeros]+1}];
Show[{ggr,Table[graphs[k],{k,Length[fpzeros]+1}]},DisplayFunction->$DisplayFunction]
]

PlotConcaveUpConcaveDown[f_,{var_,min_,max_},opts___]:=Module[{fppzeros = FindZeros[D[f,{var,2}],var],eps = 0.01 ,testvals,intervals,graphs,inc},
testvals = Join[{First[fppzeros]-eps},eps+fppzeros];
inc =Positive[D[f,{var,2}]/.{var ->testvals}];
intervals = Join[{{min,First[fppzeros]}},Table[{fppzeros[[k]],fppzeros[[k+1]]},{k,1,Length[fppzeros]-1}],{{Last[fppzeros],max}}];
Do[graphs[k] = Plot[f,{var,intervals[[k,1]],intervals[[k,2]]},opts,PlotStyle->{Thickness[0.01],If[inc[[k]],RGBColor[1,0,1],RGBColor[.4,.4,.4]]},DisplayFunction->Identity],{k,1,Length[fppzeros]+1}];
Show[Table[graphs[k],{k,Length[fppzeros]+1}],DisplayFunction->$DisplayFunction]
]

PlotConcaveUpConcaveDown[{f_,g_},{var_,min_,max_},opts___]:=Module[{fppzeros = FindZeros[D[f,{var,2}],var],eps = 0.01 ,testvals,intervals,graphs,inc,ggr},
testvals = Join[{First[fppzeros]-eps},eps+fppzeros];
inc =Positive[D[f,{var,2}]/.{var ->testvals}];
intervals = Join[{{min,First[fppzeros]}},Table[{fppzeros[[k]],fppzeros[[k+1]]},{k,1,Length[fppzeros]-1}],{{Last[fppzeros],max}}];
ggr=Plot[g,{var, min,max},PlotStyle-> {Thickness[0.01],RGBColor[0,0,1]},DisplayFunction->Identity];Do[graphs[k] = Plot[f,{var,intervals[[k,1]],intervals[[k,2]]},opts,PlotStyle->{Thickness[0.01],If[inc[[k]],RGBColor[1,0,1],RGBColor[.4,.4,.4]]},DisplayFunction->Identity],{k,1,Length[fppzeros]+1}];
Show[{ggr,Table[graphs[k],{k,Length[fppzeros]+1}]},DisplayFunction->$DisplayFunction]
]

Optimization

ShowCardboardBox[width_,length_,tocut_]:=Show[{Graphics[{{ColorData["Legacy", "Burlywood"],Rectangle[{0,0},{width,length}]},{Dashing[{0.01`,0.01`}],Red,Line[{{tocut,0},{tocut,length}}],Line[{{0,tocut},{width,tocut}}],Line[{{0,length-tocut},{width,length-tocut}}],Line[{{width-tocut,0},{width-tocut,length}}]}}],frame[{0,0},{width,length}]},AspectRatio->Automatic]
frame[corner1_,corner2_,opts__:Black]:=Graphics[{opts,Line[{corner1,{corner2[[1]],corner1[[2]]},corner2,{corner1[[1]],corner2[[2]]},corner1}]}]
CutCardboard[width_,length_,tocut_]:=Module[{gr1},gr1=Show[{Graphics[{{ColorData["Legacy", "Burlywood"],Rectangle[{0,0},{width,length}]},{Dashing[{0.01`,0.01`}],Red,Line[{{tocut,0},{tocut,length}}],Line[{{0,tocut},{width,tocut}}],Line[{{0,length-tocut},{width,length-tocut}}],Line[{{width-tocut,0},{width-tocut,length}}]}}],frame[{0,0},{width,length}]},AspectRatio->Automatic,DisplayFunction->Identity];gr2=Show[Graphics3D[Cuboid[{0,0,0},{width-2 tocut,length-2 tocut,tocut}]],AspectRatio->Automatic,Boxed->False,DisplayFunction->Identity];Show[GraphicsRow[{gr1,gr2}],DisplayFunction->$DisplayFunction,AspectRatio->Automatic,PlotLabel->"cut = "<>ToString[N[tocut]]<>"; volume = "<>ToString[N[(width-2 tocut) (length-2 tocut) tocut,8]]]]
ShowChurchWindowDetails[r_,p_]:=Module[{h=N[(1/2) (p-2 r-π r)]},Show[{Graphics[{{Pink,Disk[{r,h},r,{0,π}]},{Circle[{r,h},r,{0,π}]},{Red,Line[{{0,(h/3)},{r,(h/3)}}],Line[{{2.2` r,0},{2.2` r,h}}],Text["r",{(r/2),(h/3)},{0,1}],Text["h",{2.3` r,(h/2)},{-1,0}],Text["r",{r+(((Sqrt[2.`]) r)/4),h+(((Sqrt[2]) r)/4)},{0,1}],Line[{{r,h},{r+(((Sqrt[2.`]) r)/2),h+(((Sqrt[2]) r)/2)}}]}}],frame[{0,0},{2 r,h}]},AspectRatio->Automatic]]
ShowChurchWindow[r_,p_]:=Module[{h=(1/2) (p-2 r-π r)},Show[{Graphics[{{Pink,Disk[{r,h},r,{0,π}]},{Circle[{r,h},r,{0,π}]}}],frame[{0,0},{2 r,h}]},AspectRatio->Automatic]]
ShowFence[width_,fencing_]:=Module[{h=((fencing-width)/2)},Show[Graphics[{{YellowGreen,Rectangle[{-(fencing/2),0},{(fencing/2),fencing}]},{AbsoluteThickness[2],Line[{{-(width/2),h},{(width/2),h}}],Line[{{-(width/2),0},{-(width/2),h}}],Line[{{(width/2),0},{(width/2),h}}]},{Blue,Rectangle[{-(fencing/2),-20},{(fencing/2),0}]},{Green,Text["River",{0,-10}]}}],PlotRange->All,AspectRatio->Automatic,PlotLabel->"width = "<>ToString[width]<>"\narea = "<>ToString[N[width h]]]];
ShowRectangle[w_,r_]:=Module[{h=N[(Sqrt[(r^2)-((w^2)/4)])]},Show[Graphics[{{Blue,Circle[{0,0},r,{0,π}],Line[{{-r,0},{r,0}}]},{Red,Rectangle[{-(w/2),0},{(w/2),h}]}},AspectRatio->Automatic,PlotLabel->"width = "<>ToString[w]<>"\narea = "<>ToString[w h]]]]
ShowCuts[fn_,{var_,min_,max_},n_Integer,opts___]:=ShowCuts[fn,{var,min,max},Table[i,{i,min,max,((max-min)/n)}],opts]
ShowCuts[fn_,{var_,min_,max_},pts_List,opts___]:=Module[{xpt,fnx,gr1,gr2},fnx=Function[var,fn];gr1=Show[Graphics[{Blue,Table[xpt=pts[[i]];Line[{{xpt,fnx[xpt]},{xpt,-fnx[xpt]}}],{i,Length[pts]}]}],DisplayFunction->Identity,AspectRatio->Automatic];gr2=FilledPlot[{fn,-fn},{var,min,max},Fills->Yellow,PlotStyle->Black,AspectRatio->Automatic,DisplayFunction->Identity,opts];Show[{gr2,gr1},AspectRatio->Automatic,DisplayFunction->$DisplayFunction,opts]]

ParametricPlot

TraceParametricPlot[{ff_,gg_},{t_,t0_,t1_},dt_Integer:8,opts___]:=TraceParametricPlot[{f,g},{t,t0,t1,dt},opts]

TraceParametricPlot[{ff_,gg_},{t_,t0_,t1_,dt_Integer:8},opts___]:=Module[{ndt,arro,arrsz,rng,xx,tt,f,g},ndt=N[(t1-t0)/dt];
f=Function[t,ff];g=Function[t,gg];
arrsz=N[(Max[Table[Abs[f[tt]],{tt,t0,t1,N[(t1-t0)/25]}]]+Max[Table[Abs[g[tt]],{tt,t0,t1,N[(t1-t0)/25]}]])/2.0];
arro[{a_,b_},{c_,d_},eps_]:=Polygon[{{a,b}+3*eps*{c-a,d-b}/Sqrt[(b-d)^2+(c-a)^2],{a,b}+eps/Sqrt[(b-d)^2+(c-a)^2]*{(b-d),(c-a)},{a,b}+(-eps)/Sqrt[(b-d)^2+(c-a)^2]*{(b-d),(c-a)}}];
rng=FullOptions[ParametricPlot[{f[xx],g[xx]},{xx,t0,t1},opts,PlotRange->All,DisplayFunction->Identity],PlotRange];
Do[ParametricPlot[{f[xx],g[xx]},{xx,t0,tt},PlotRange->rng+{{0,0},{0,.06*rng[[2,2]]}},PlotStyle->Hue[(tt+1)/(dt-1)],Epilog->{Blue,arro[{f[tt],g[tt]},{f[tt+.05*arrsz],g[(tt+.05*arrsz)]},.035*arrsz]},PlotLabel->StringJoin["t = ",ToString[N[If[Abs[tt+.001]<.0001,0,tt+.001],3]]],DisplayFunction->$DisplayFunction,opts],{tt,t0-.001,t1,ndt}]]

RainbowParametricPlot[{fx_,fy_},{var_,min_,max_},n_Integer:8,opts___]:=Module[{dx=(max-min)/n,gr},
Do[gr[j]=ParametricPlot[{fx,fy},{var,min+(j-1) dx,min+(j) dx},PlotStyle->Hue[j/n],opts],{j,1,n}];
Show[Table[gr[i],{i,1,n}], PlotRange->All, AxesOrigin->{0,0}]]

RainbowParametricPlotwDots[{ff_,gg_},{var_,min_,max_},n_Integer:8,opts___]:=Module[{ln,gr,rt,i,pt,rng,dx,st,grs,f,g,pts={}},dx=(max-min)/n;
f=Function[var,ff];
g=Function[var,gg];
Do[st=min+(i-1) dx;
AppendTo[pts,{PointSize[0.01+i 0.0025],Point[{f[st],g[st]}]}];
gr[i]=ParametricPlot[{f[var],g[var]},{var,st,st+dx},PlotStyle->Hue[i/n],opts],{i,1,n}];grs=Table[gr[i],{i,1,n}];
Show[{grs,Graphics[pts]}, PlotRange->All, AxesOrigin->{0,0}]]

RainbowPolarPlot[r_,{var_,min_,max_},n_Integer,opts___]:=Module[{ln,gr,rt,i,pt,rng,dx,st,grs},dx=(max-min)/n;
rt=Function[var,r];
Do[st=min+(i-1)dx;
gr[i]=PolarPlot[r,{var,st,st+dx},PlotStyle->Hue[i/n]],{i,1,n}];
grs=Table[gr[i],{i,1,n}];
Show[grs]]

TracePolarPlot[r_,{var_,min_,max_},n_Integer,opts___]:=TracePolarPlot[r,{var,min,max},Table[i,{i,min+(max-min)/n,max,(max-min)/n}],opts]

TracePolarPlot[r_,{var_,min_,max_},pts_List,opts___]:=Module[{n,ln,gr,rt,i,pt,rng},n=Length[pts];
rt=Function[var,r];
rng=1.2*FullOptions[PolarPlot[r,{var,min,max},opts,DisplayFunction->Identity],PlotRange];
Do[pt=pts[[i]];
gr=PolarPlot[r,{var,min,pt},DisplayFunction->Identity];
ln=Graphics[{Blue,Line[{{0,0},{trigcossign[rt[pt] Cos[pt],pt],trigsinsign[rt[pt] Sin[pt],pt]}}]}];
Show[{gr,ln},DisplayFunction->$DisplayFunction,PlotLabel->StringJoin[ToString[N[pt,5]]],PlotRange->rng,opts],{i,1,n}]]

Series

7.5 Approximating Definite Integrals

Consider the function f(x)=(Sqrt[x]) and the area underneath it from x=1 to x=5. Let's divide the interval into 4 parts (for now, using num=4) and look at approximations to the area using left, right, and midpiont-based rectangles, as well as a trapezoid-based approach.

num = 4;
f[x_]:=(Sqrt[x])

ShowApproximations[f[x],{x,1,5},left,num,PlotRange->{{0,5},Automatic}]

ShowApproximations[f[x],{x,1,5},right,num,PlotRange->{{0,5},Automatic}]

ShowApproximations[f[x],{x,1,5},midpoint,num,PlotRange->{{0,5},Automatic}]

ShowApproximations[f[x],{x,1,5},trapezoid,num,PlotRange->{{0,5},Automatic}]

Note that the function above is an increasing function. Did the left-based approach underestimate or overestimate the actual value? Could you tell how things were with the midpoint-based approach? How about the trapezoidal approach?
Now we do the same thing again but with a function that is still monotone, but decreasing.

f[x_]:=(1/x)

ShowApproximations[f[x],{x,1,5},left,num,PlotRange->{{0,5},Automatic}]

ShowApproximations[f[x],{x,1,5},right,num,PlotRange->{{0,5},Automatic}]

ShowApproximations[f[x],{x,1,5},midpoint,num,PlotRange->{{0,5},Automatic}]

ShowApproximations[f[x],{x,1,5},trapezoid,num,PlotRange->{{0,5},Automatic}]

Answer the same questions posed above.

Now consider the following four functions. Think about increasing/decreasing and concavity and classify these.

f[x_]:= (x^2)
g[x_]:=(Sqrt[x])
h[x_]:= (1/x)
k[x_]:=-(x^2)+36

Without seeing the graphs, can you classify each of the preceding as increasing or decreasing and as concave up or concave down. What implications do these attributes, if any, have on whether we have an overestimate or underestimate? Let's investigate just the trapezoidal approach for now.

num = 2;
gr1=ShowApproximations[f[x],{x,1,5},trapezoid,num,PlotRange->{{0,5},Automatic}]
gr2=ShowApproximations[g[x],{x,1,5},trapezoid,num,PlotRange->{{0,5},Automatic}]
gr3=ShowApproximations[h[x],{x,1,5},trapezoid,num,PlotRange->{{0,5},Automatic}]
gr4=ShowApproximations[k[x],{x,1,5},trapezoid,num,PlotRange->{{0,5},Automatic}]

Here, we put them all together.

Show[GraphicsGrid[{{gr1,gr2},{gr3,gr4}}]]

Now we investigate the midpoint approach.

num = 2;
gr1=ShowApproximations[f[x],{x,1,5},midpoint,num,PlotRange->{{0,5},Automatic}]
gr2=ShowApproximations[g[x],{x,1,5},midpoint,num,PlotRange->{{0,5},Automatic}]
gr3=ShowApproximations[h[x],{x,1,5},midpoint,num,PlotRange->{{0,5},Automatic}]
gr4=ShowApproximations[k[x],{x,1,5},midpoint,num,PlotRange->{{0,5},Automatic}]

Show[GraphicsGrid[{{gr1,gr2},{gr3,gr4}}]]

These may be a little difficult to see. However, every midpoint approximation is really a trapezoidal one using the tangent line to the curve at the midpoint. Study what is shown below and see if you can see this. Also, answer the question that we have been asking about how increasing/decreasing and/or concavity effect being an overestimate or underestimate.

num = 2;
gr1=ShowApproximations[f[x],{x,1,5},modmidpoint,num,PlotRange->{{0,5},Automatic}]
gr2=ShowApproximations[g[x],{x,1,5},modmidpoint,num,PlotRange->{{0,5},Automatic}]
gr3=ShowApproximations[h[x],{x,1,5},modmidpoint,num,PlotRange->{{0,5},Automatic}]
gr4=ShowApproximations[k[x],{x,1,5},modmidpoint,num,PlotRange->{{0,5},Automatic}]

Show[GraphicsGrid[{{gr1,gr2},{gr3,gr4}}]]

7.6 Approximation Errors and Simpson’s Rule

A quick review

Consider the function f(x) = -2(x+1)(x-1)(x-3). Suppose we are interested in the area of the region underneath the curve, above the x-axis and between x = 1 and x = 3, as shown below:

Clear[f,x]
f[x_]:=-2 (x+1) (x-1) (x-3)
Plot[f[x],{x,1,3},Filling->Axis, FillingStyle ->RGBColor[1,1,0]]

Earlier we have considered various rectangular approximations, as well as trapezoidal approximations. The rectangular types include: left, right, midpoint, random, inscribed, and circumscribed. Below, if you wish, you can review how these rectangular approximations worked.

ShowApproximations[f[x],{x,1,3},left,6]

ShowApproximations[f[x],{x,1,3},right,6]

ShowApproximations[f[x],{x,1,3},midpoint,6]

ShowApproximations[f[x],{x,1,3},random,6]

ShowApproximations[f[x],{x,1,3},inscribed,6]

ShowApproximations[f[x],{x,1,3},circumscribed,6]

Here, we will focus on the first three type of rectangular approximations, the trapezoidal approximation and the Simpson method. We will learn more about the Simpson Rule later in this lab. When we wish to refer to an approximation using n subdivisions (of the domain), we will use the notation (Subscript[L, n]), (Subscript[R, n]), (Subscript[M, n]), (Subscript[T, n]), and (Subscript[S, n]) to represent the estimates of the left, right, midpoint, trapezoid and Simpson methods, respectively. The cell below yields numerical estimates using these five methods.

LeftApprox[f[x],{x,1,3},8]
RightApprox[f[x],{x,1,3},8]
MidpointApprox[f[x],{x,1,3},8]
TrapezoidApprox[f[x],{x,1,3},8]
SimpsonApprox[f[x],{x,1,3},8]

To see the results in decimal form, use the cousins of each of these with an initial N appended:

NLeftApprox[f[x],{x,1,3},8]
NRightApprox[f[x],{x,1,3},8]
NMidpointApprox[f[x],{x,1,3},8]
NTrapezoidApprox[f[x],{x,1,3},8]
NSimpsonApprox[f[x],{x,1,3},8]

(For the ten commands just illustrated, each one has an alternate name with the ‘Approx’ dropped and all lowercase letters used.) The actual value can be obtained by integrating.

Integrate[f[x], {x, 1, 3}]

In this case, we see that the Simpson method yielded an exact value.

Q1. What was the next best method (in this case)?

Answer:

Let’s consider another function, f(x) = sin(2x). Suppose we are interested in the area of the region underneath the curve, above the x-axis and between x = 0 and x = 5π/12, as shown below:

Clear[f,x]
f[x_]:=Sin[2 x]
Plot[f[x],{x,0,((5 π)/12)}]

Evaluate the cell below to obtain estimates.

LeftApprox[f[x],{x,0,((5 π)/12)},8]
RightApprox[f[x],{x,0,((5 π)/12)},8]
MidpointApprox[f[x],{x,0,((5 π)/12)},8]
TrapezoidApprox[f[x],{x,0,((5 π)/12)},8]
SimpsonApprox[f[x],{x,0,((5 π)/12)},8]

It is clear that in this case, the decimal versions will be more useful.

NLeftApprox[f[x],{x,0,((5 π)/12)},8]
NRightApprox[f[x],{x,0,((5 π)/12)},8]
NMidpointApprox[f[x],{x,0,((5 π)/12)},8]
NTrapezoidApprox[f[x],{x,0,((5 π)/12)},8]
NSimpsonApprox[f[x],{x,0,((5 π)/12)},8]

Q2. Determine the exact value by integrating (either by hand or by using Mathematica). Now order the five results above from best to worst. Explain how you determined your ranking.

Answer:

N[(∫_0^((5 π)/12))f[x]\[DifferentialD]x]

Error analysis of monotone functions

Consider the following two graphs of the function h(x) = (x^2) on the interval [1, 3].

Clear[h,x]
h[x_]:=(x^2)
ShowApproximations[h[x],{x,1,3},left,8]
ShowApproximations[h[x],{x,1,3},right,8]

Here we see that (Subscript[L, 8])=7.6875 while (Subscript[R, 8])=9.6875. An important observation can be made here. Since the function h is monotone (in this case increasing) on the interval, then the left endpoint method actually results in inscribed rectangles. This means the estimate is guaranteed to be a lower bound on the actual value of the integral. Similarly, the right endpoint method results in circumscribed rectangles, resulting in an upper bound of the actual value. Let’s consider this more closely.

ShowLeftRightDifference[h[x],{x,1,3},8]

When you look at the graph above, what you see in blue represents the difference between the left and right estimates.

Q3. What is the area of these blue rectangles? Note that they each have the same width. Furthermore, observe that the height of the top of the ith rectangle is the height of the bottom of the (i+1)st rectangle. Therefore, imagine sliding them all together on one stack (say, on the first rectangle). How great is this height? Consequently, what is the area?

Answer:

How much error is there with the trapezoidal method? Consider the following.

ShowApproximations[h[x],{x,1,3},trapezoid,8]

This method is quite good. If we reduce the number of subintervals and/or change functions, the error is more apparent. The second graph uses red to show the “missed” regions.

ShowApproximations[Sin[x],{x,1,3},trapezoid,2]
ShowTrapError[Sin[x],{x,1,3},2]

Consider the following theorem (as found in the Ostebee/Zorn text, Calculus from Graphical, Numerical and Symbolic Points of View, section 7.1). Recall from the previous section the definitions of (Subscript[R, n]), (Subscript[L, n]), and (Subscript[T, n]).

Theorem 1: Suppose that f is monotone on [a, b] and I =  (∫_a^b)f(x)dx. Then:

|I-(Subscript[R, n])|≤|f(b)-f(a)|((b-a)/n)

|I-(Subscript[L, n])|≤|f(b)-f(a)|((b-a)/n)

|I-(Subscript[T, n])|≤|f(b)-f(a)|((b-a)/(2n))

Let’s consider this theorem with our function h(x) = (x^2) on [1, 3]. First let’s calculate what the integral is exactly. This is easy to do by hand (or in your head), but we will use Mathematica.

Integrate[(x^2),{x,1,3}]

We have already determined (Subscript[L, 8]), (Subscript[R, 8]), and (Subscript[T, 8]), so we can calculate the actual errors.

lefterror=Abs[(26/3)-7.6875]
righterror=Abs[(26/3)-9.6875]
trapezoiderror=Abs[(26/3)-8.6875]

Are these actual errors within the bounds guaranteed by the theorem? (They had better be, or our theorem is not really a theorem!) Let’s compute these bounds.

n=8;
a=1;
b=3;
((Abs[h[b]-h[a]] (b-a))/n)
((Abs[h[b]-h[a]] (b-a))/(2 n))

The theorem guarantees the error to be less than 2 for left or right endpoints and less than 1 for the trapezoid method; our results indeed fall within these bounds.

Q4. Using the same function, increase the number of subintervals to 16. What are the actual errors and what are the error bounds?

Answer:

Q5. Now consider the function k(x) = sin(x) on the interval [0, π] and suppose the number of subintervals is 10. What are the actual errors and what are the error bounds in this case? Does the theorem apply? Why or why not?

Answer:

Error analysis, in general, for left and right methods

Consider the constant function p(x) = 4 on the interval [1, 5]. What will the left endpoint method yield for an approximation and what is the exact area. Think about this before proceeding. OK, now let’s use our tools to confirm our thoughts.

Clear[p,x]
p[x_]:=4
ShowApproximations[p[x],{x,1,5},left,6]

Hopefully you thought about the result before you evaluated the above cell and wondered why we even should ask for help from Mathematica. The picture should make it blatantly clear that for constant functions, the left (and right) method has no error. Now consider the following two examples:

ShowApproximations[(x/2),{x,1,5},left,6,PlotRange->{(1/2),10},AxesOrigin->{1,0}]
ShowApproximations[2 x,{x,1,5},left,6,PlotRange->{(1/2),10},AxesOrigin->{1,0}]

Note that the y-axes are scaled the same on both graphs.

Q6. Which function has more error? Why?

Answer:

How would the error of the function y = 10 x compare to these last two? Note that the only difference in the two examples above is the slope of the function. The slope, of course, is measured by the first derivative. Let’s consider a table of numerical data.

functlist={(x/10)-4,(x/4)+2,x,5 x-3,-9 x+3,15 x,-23 x+1}
TableForm[Table[{functlist[[i]],est=NLeftApprox[functlist[[i]],{x,0,4},6],act=N[(∫_0^4)functlist[[i]]\[DifferentialD]x],error=act-est,Abs[error]},{i,1,Length[functlist]}],TableHeadings->{None,{"f","L(6)","I = exact","I - L(6)","|I - L(6)|"}},TableAlignments->{Center,Center}]

The first column is a list of some linear functions. The second column is (Subscript[L, 6]) for each function over the interval [0, 4]. The third column is the exact integral and the next two represent the error and absolute error. Think about this table and observe what happens to the error as we change the slope of the line. Consider also the following table which considers the function y = x on the interval [0, 4] and considers increasing the number of subintervals (n).

TableForm[Table[{(2^i),est=NLeftApprox[x,{x,0,4},(2^i)],act=N[(∫_0^4)x\[DifferentialD]x],Abs[act-est]},{i,0,6}],TableHeadings->{None,{"n","L(n)","I = exact","|I - L(n)|"}},TableAlignments->{Center,Center}]

Q7. Examine the two tables closely. What conclusions can you make regarding the error for the left endpoint method?

Answer:

A theorem

Consider the following theorem (as found in the Ostebee/Zorn text, Calculus from Graphical, Numerical and Symbolic Points of View,, section 7.2).

Theorem 2: Let I = (∫_a^b)f(x)dx and let (Subscript[L, n]) be the left approximating sum for I, with n equal subdivisions. Suppose that for all x in [a, b], |f'(x)|(Subscript[≤K, 1]). Then

|I-(Subscript[L, n])|≤((((Subscript[K, 1])(b-a))^2)/(2n)).

(The same error holds for right sums.)

To understand this theorem, you need to know more about (Subscript[K, 1]). The questions above should have pointed you in the direction of seeing a relationship between the amount of error in the left approximation and the slope of the function (lines in the cases that we examined). We use (Subscript[K, 1]) as a value of an upper bound for |f'| (the absolute value of the first derivative---recall that the derivative measures slope) on the interval [a, b]. Let’s consider an example. Suppose we consider f(x) = -2(x+1)(x-1)(x-3), which was our first example in the first section.

Clear[f,x,n,a,b]
f[x_]:=-2 (x+1) (x-1) (x-3)
n=6;a=1;b=3;
ShowApproximations[f[x],{x,a,b},left,n]

In this case, there is internal and external error (ie., some regions under the curve were missed and some regions outside the region were covered). Now consider f' and its graph:

f'[x]
Simplify[%]

Plot[f'[x],{x,a,b}]

From this, it appears as if |f'(x)| < 18 for all x in the interval [1, 3]. So, we could use (Subscript[K, 1])=18. (Perhaps we could use 16.5, but 18 is also an upper bound.) Note that the graph of |f'| might be slightly easier to read.

Plot[Abs[f'[x]],{x,a,b}]

Using (Subscript[K, 1])=18, from the theorem we can calculate the upper bound for our error:

K1=18;
maxerror=((K1 (((b-a))^2))/(2 n))

Q8. We just determined that the maximum the error could be is 6. What was the actual error?

Answer:

two more examples

Let’s consider another function. Often one hears the phrase “graded on the curve” or something similar. What curve? Well, consider f(x) = (1/(Sqrt[2π]))(e^(((-(x^2))/2) )) on the interval [-4, 4] with 8 subdivisions:

Clear[f,x,a,b,n]
f[x_]:=((E^(-((x^2)/2)))/(Sqrt[2 π]))
a=-4;b=4;n=8;
Plot[f[x],{x,a,b}]
ShowApproximations[f[x],{x,a,b},left,n]

The left estimate yields about 0.999863. How far away is this from the exact value? Let’s use our theorem, which gives the bounds for this error. First we need to determine (Subscript[K, 1]), so let’s plot |f'|.

Plot[Abs[f'[x]],{x,-4,4}]

It looks like we could safely say (Subscript[K, 1])=0.25. Now to calculate the error bound.

K1=0.25;
maxerror=((K1 (((b-a))^2))/(2 n))

Since we obtained 0.999863 for (Subscript[L, 8]), this means that |I - 0.999863| ≤ 1, or we have -1 ≤ I - 0.999863 ≤ 1, or -0.000137 ≤ I ≤ 1.999863. This doesn’t seem like much help in obtaining an estimate of I. Let’s increase our value of n to 100.

n=100;
NLeftApprox[f[x],{x,a,b},n]

Thus, we see that (Subscript[L, 100])=0.999936. We should be fairly confident of this, since we used so many intervals. Let’s determine the error bounds. Note that (Subscript[K, 1]) doesn’t change, only n.

K1=0.25;
maxerror=((K1 (((b-a))^2))/(2 n))

So now we can be sure that the exact value is 0.999936 ± 0.08. So let’s calculate the exact value:

Integrate[f[x], {x, a, b}]

What in the world is this? This does NOT look a numerical value. Can’t Mathematica compute the antiderivative of our function and apply the Fundamental Theorem of Calculus? The truth is that not only it can’t, but neither can you---there is no elementary antiderivative! (Actually, Erf does mean something, but we won’t worry about it.) We could ask Mathematica to use its numerical approach to integration by using NIntegrate.

NIntegrate[f[x], {x, a, b}]

Of course, this is also done by estimation.

Let’s consider another function which does not have an (elementary) antiderivative: k(x) = cos((x^2)) (and use right sums instead of left—the theorem applies to both)

Clear[k,a,b,x,n]
k[x_]:=Cos[(x^2)]
a=-1;b=1;n=20;
Plot[k[x],{x,a,b}]
ShowApproximations[k[x],{x,a,b},right,n]
NRightApprox[k[x],{x,a,b},n]

So an estimate is 1.80624. What will Integrate yield?

(∫_a^b)k[x]\[DifferentialD]x

Although Mathematica does return something when we ask it to integrate our function, it is not a numerical value of something we recognize.

Q9. Determine the error bounds for this function. What do you think is the “actual” value? How did you arrive at this estimate?

Answer:

Above we used right sum estimates for k(x) = cos((x^2)) on [-1, 1] with 20 subdivisions. Consider the other possibilities:

NRightApprox[k[x],{x,a,b},n]
NLeftApprox[k[x],{x,a,b},n]
NMidpointApprox[k[x],{x,a,b},n]
NTrapezoidApprox[k[x],{x,a,b},n]

Q10. Do these results effect what you think might be the “actual” value?

Answer:

Now your turn

Consider the function f(x) = (Sqrt[x]) on the interval [1, 3]. The theorem for error bounds gives us an upper bound on the error using n intervals. We know that this involves a, b, n, and (Subscript[K, 1]). Suppose we wanted a value which is accurate to 1 decimal place. In other words, we must have |I-(Subscript[L, n])|<0.05. But then we must have 
((((Subscript[K, 1])(b-a))^2)/(2n))<0.05 to ensure |I-(Subscript[L, n])|<0.05. What value of n will do this? First, let’s determine (Subscript[K, 1]):

Clear[f,x,a,b,n]
f[x_]:=(Sqrt[x])
a=1;b=3;
Plot[Abs[f'[x]],{x,a,b},PlotRange->All]

It appears as if we could safely use (Subscript[K, 1])=0.55. (Remember, this only has to be some upper bound, but in this case we would like a bound fairly close to the least upper bound.) Let’s try this out.

K1=0.55
maxerror=((K1 (((b-a))^2))/(2 n))

Since we haven’t specified n (and we included it in our Clear statement to prevent carryover from a previous problem), we must expect this to be a function of n. We want this to be less than 0.05. Of course, this is simple enough to calculate by hand, but let’s see how we can approach this graphically. In other words, when does this function of n intercept our constant 0.05? Certainly n > 1.

Plot[{0.05`,maxerror},{n,1,20}]

Not far enough out; try again.

Plot[{0.05`,maxerror},{n,21,40}]

Approximately, n=23. What does this mean? If we use at least 23 subdivisions, we can be guaranteed that our error will be such that we will have one decimal of accuracy. Is this true in practice? The exact value can be easily calculated.

exact = Integrate[f[x], {x, a, b}]//N

Since we tolerate error up to 0.05, we can determine the interval in which we will be satisfied that we have one decimal accuracy.

tolerance=0.05;
interval={exact-tolerance,exact+tolerance}

So, if our estimate lands in this interval of “satisfaction,” we will be content that the value has accuracy to the first decimal. Let’s make a table to calculate some left approximations:

TableForm[Table[{i,NLeftApprox[f[x],{x,a,b},i]},{i,5,15}]]

From this, we see that we enter the interval [2.74743, 2.84743] when we use n = 15. In summary, the theorem indicates that if we use 23 subintervals we can be guaranteed to be within 0.05 of the actual value, but in practice we only need 15 subintervals.

Q11. Suppose we wish to know the area under the sine curve on the interval [0, π] with one decimal accuracy (ie., be within 0.05 of the actual value). Modify the steps above to determine the number of subintervals the theorem guarantees will accomplish this and also determine the smallest value of n which is actually needed.

Answer:

Error analysis for midpoint and trapezoid methods

We now wish to consider the midpoint and trapezoid methods. First let’s look at examples of both:

Clear[h,x,a,b,n]
h[x_]:=(x^2)
a=-3;b=3;n=8;
ShowApproximations[h[x],{x,a,b},trapezoid,n]
ShowApproximations[h[x],{x,a,b},midpoint,n]
{"M(30) = "<>ToString[NMidpointApprox[h[x],{x,a,b},30]],"T(30) = "<>ToString[NTrapezoidApprox[h[x],{x,a,b},30]]}

This provides some evidence that these two are quite close to each other. The graphic below is a “Proof Without Words” illustrating that the midpoint method can be considered very similar to the trapezoid method.

RectToTrap

Q12. Explain what these pictures are trying to show.

Answer:

Whereas with the left and right methods, the slope  of the curve effected the amount of error, with these two methods the concavity effects the error. Accordingly, define (Subscript[K, 2]) to be an upper bound of |f''| on the interval [a, b]. We then have the following theorem (as found in the Ostebee/Zorn text, Calculus from Graphical, Numerical and Symbolic Points of View,, section 7.3):

Theorem 3: Let I = (∫_a^b)f(x)dx and let (Subscript[K, 2]) to be an upper bound of |f''| on the interval [a, b]. Then

|I-(Subscript[M, n])|≤ (( (((Subscript[K, 2])(b-a))^3))/(24(n^2))) and |I-(Subscript[T, n])|≤ (( (((Subscript[K, 2])(b-a))^3))/(12(n^2))).

Let’s apply this to a familiar function: f(x) = -2(x+1)(x-1)(x-3) on [1, 3]

Clear[f,x,a,b,n]
f[x_]:=-2 (x+1) (x-1) (x-3)
a=1;b=3;n=8;
ShowApproximations[f[x],{x,a,b},trapezoid,n]
ShowApproximations[f[x],{x,a,b},midpoint,n]
(∫_a^b)f[x]\[DifferentialD]x

We see that the exact value is 8, while the midpoint method overestimated with 8.0625 and the trapezoid method underestimated with 7.875, yielding errors of 0.0625 and 0.125 respectively. What does the theorem guarantee for the largest amount of error that can occur when n = 8?

Plot[Abs[f''[x]],{x,a,b}]

It appears as if we could use (Subscript[K, 2])=25. Therefore, our errors are guaranteed to be less than the following:

K2=25;
miderror=N[((K2 (((b-a))^3))/(24 (n^2)))]
traperror=N[((K2 (((b-a))^3))/(12 (n^2)))]

Note that both of these are significantly higher than what we actually experienced.

Q13. Using the same function and interval, what are the bounds for errors when 30 subintervals are used?

Answer:

Q14. This question is similar to number 11 in the previous section. Suppose we wish to know the area shown below (evaluate the cells) with one decimal accuracy (ie., be within 0.05 of the actual value). Determine the number of subintervals the theorem guarantees will accomplish this (using either the trapezoid or midpoint method) and also determine the smallest value of n which is actually needed.

Answer:

Clear[f,x,a,b,n]
f[x_]:=895-1404 x+1282 (x^2)-976 (x^3)+499 (x^4)-139 (x^5)+19 (x^6)-(x^7)
a=1;b=6;
Plot[f[x],{x,1,6}]

The Simpson Rule

We now wish to explore another method of approximation---Simpson’s rule. Consider the same function used in the last exercise:

Clear[f,x,a,b,n]
f[x_]:=895-1404 x+1282 (x^2)-976 (x^3)+499 (x^4)-139 (x^5)+19 (x^6)-(x^7)
a=1;b=6;n=10;
ShowApproximations[f[x],{x,a,b},simpson,n]

What does this graph show? With Simpson’s rule, we MUST have an EVEN number of subintervals. Here, we used 10 subdivisions. Recall that the left and right rectangle methods used a degree 0 polynomial (constant function) for a “top,” while the midpoint and trapezoid methods used a degree 1 polynomial (linear function) as a top. With Simpson’s rule we use a degree 2 polynomial (quadratic function) top. Here we subdivide the interval into an even number of subintervals. For each successive pair of subintervals, we determine the unique parabola which goes through the three points determined by the pair of subintervals. (Look at the graph above---a dotted line corresponds to the midpoint of the pair of subintervals.) We then draw the parabolas in red and the original function is drawn in black. Note the error on the interval [5, 6].

How good is the Simpson method? Let’s use the same function as above and compare various methods with several values of n.

TableForm[Table[{n,NLeftApprox[f[x],{x,a,b},n],NRightApprox[f[x],{x,a,b},n],NTrapezoidApprox[f[x],{x,a,b},n],NMidpointApprox[f[x],{x,a,b},n],NSimpsonApprox[f[x],{x,a,b},n]},{n,8,30,4}],TableHeadings->{None,{"n","L(n)","R(n)","T(n)","M(n)","S(n)"}},TableAlignments->{Center,Center}]

With the other methods, we saw that certain bounds on the derivative effected the error bounds. In this case, define (Subscript[K, 4]) (you might ask “What happened to (Subscript[K, 3])?”) to be an upper bound of |f''''| on the interval [a, b]. We then have the following theorem (as found in the Ostebee/Zorn text, Calculus from Graphical, Numerical and Symbolic Points of View,, section 7.4):

Theorem 4: Let I = (∫_a^b)f(x)dx and let (Subscript[K, 4]) to be an upper bound of |f''''| on the interval [a, b]. Then

|I-(Subscript[S, n])|≤((((Subscript[K, 4])(b-a))^5)/(180(n^4))).

Q15. Repeat question 13 using the Simpson rule theorem.

Answer:

Q16. Repeat question 14 using the Simpson rule theorem.

Answer:

1. Improper Integral example: integrable function with a vertical asymptote

One of our questions from this week’s material asks “What makes (∫_0^1)(1/(Sqrt[x]))d x an improper integral? Does it converge or not?”
It’s improper because of the vertical asymptote at the left endpoint.
It does converge, because it blows up “slowly enough” that there is actually a finite amount of area under the curve, even as the curve itself heads to infinity.

Define f(x)=(1/(Sqrt[x])), a and b as the endpoints.

Use the right endpoint method to approximate this integral with n=100, using Mathematica’s Sum command. Here it a picture of the situation with n=8: :

Use the midpoint method to approximate this integral with n=100, using Mathematica’s Sum command. Here it a picture of the situation with n=8:

Why are you unable to use any other of our methods for this definite integral?

Why are you unable to carry out any error analysis using our error formula for midpoint?

Nonetheless, you can compute this integral analytically, and so compare the absolute error of both the midpoint and right endpoint methods.

2. Numerical Integration with data:

An Exercise in Modeling Human Temperature

The Data

The following figure is reproduced from the Book “La Chaleur Animale” (Animal Heat), 1889, by Charles Richet (Nobel Prize in Medicine, 1913). 


I digitized the data (pulled it off that figure) as follows (and you might check me!). The following commands define the data:

temps=37+ (* degrees celsius *)
.1{0,-1,-2,-3,-3,-3,-3,-2,-2,0,2,3,3,4,4,4,5,5,5,4,4,3,2,1}
times={0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5,11.5,12.5,13.5,14.5,15.5,16.5,17.5,18.5,19.5,20.5,21.5,22.5,23.5}(* half hours from midnight *)
data=Transpose[{times,temps}];

The model

A model, obtained by non-linear regression is
37.125` -0.40418159351544297` Sin[(1/12) π (1.0635459627341117` +t)]
and a plot of the data and the model will look like this: 


1. Define nlm(t) by the formula above (nlm stands for non-linear model).

If you defined nlm(t) properly, then you’ll generate the graph above with this command:

Show[ListPlot[data],Plot[nlm[t],{t,0,24}],PlotLabel-> "Human Temperature Across a Day",Frame->True]

The mean obtained by the non-linear regression model is 37.125 (in degrees Celsius). 
	a. What is this temperature in degrees Fahrenheit? 
	b. At what time of day is the temperature highest, according to the model? 
	c. At what time does the body temperature equal its mean?
	d. How much does the temperature vary across the day (in degrees C)?

2. Use the original data and the Trapezoidal Rule to compute the average value of temperature. How does it compare to the average obtained by the model?

3. An example problem from statistics:

This function cannot be evaluated analytically by any of our methods -- the famous “bell curve”, so essential for statistics:

f[x_]:=E^(-(x^2))
Plot[f[x],{x,-2,2}]

Our objective is to compute the integral A=(∫_(-1)^1)f(x) d x.

Compute the integral  (∫_0^1)f(x) d x using Mathematica, and call this value true (Even though Mathematica is also doing a numerical approximation!).

Why can we focus all of our energy and attention on the interval [0,1]?

Graph the absolute values of the second and fourth derivatives, and obtain estimates for the error constants K for midpoint, trapezoidal and Simpson’s rules.

Given the following estimates for midpoint ((Subscript[M, 8])), trapezoidal ((Subscript[T, 8])), and Simpson’s ((Subscript[S, 16])) rules.

NMidpointApprox[f[x],{x,0,1},8]
NTrapezoidApprox[f[x],{x,0,1},8]
NSimpsonApprox[f[x],{x,0,1},16]

Verify that the actual absolute errors of each of the three methods satisfy the error bounds associated with the given numbers of subintervals on [0,1].

Verify that Simpson’s approximation is the appropriate average of the midpoint and trapezoidal approximations.

How do we extend both the approximations and the error estimates to A, over the entire interval [-1,1]?

4. Error Analysis: you compute the number of subdivisions needed.

Compute the integral  (∫_0^4)(x^3)  (Sqrt[(x^2)+9 ])d x

analytically.

Determine what the minimal value of n should be for the midpoint rule to have accuracy within 0.0001, and compute the midpoint estimate with that value of n. Verify that the absolute error of your estimate is within the error bound.

Determine what the minimal value of n should be for Simpson’s rule to have accuracy within 0.0001, and compute the Simpson’s rule  estimate with that value of n. Verify that the absolute error of your estimate is within the error bound.