leftbox.mws
Animations for the Calculus Classroom:
Definite Integral Definitions, Approximations and Applications
by Alan Gold, Mathematics & Statistics, University of Windsor, Canada e-mail: gold1@uwindsor.ca
Note: This worksheet provides a number of plot procedures and animations for instructors to use in explaining the Definite Integral, its approximations and aplications. They are not intended as procedures for students to learn themselves, except perhaps as "black boxes".
Problem and Objectives: The Maple Student library provides graphical illustrations of approximations to definite integrals using the
leftbox
,
rightbox
and
middlebox
commands. This worksheet first shows (Section A) how to incorporate these commands into an animation to illustrate convergence of sums of rectangular areas to areas under a curve, thus illustrating the Riemann integral for these particular cases. The worksheet next provides additional procedures for the more general situations where the sample point is randomly chosen in each cell of a uniform partition (Section B), and finally where the partitions themselves are randomly generated subject to a bound on the length of the longest subinterval (Section C). The worksheet continues with procedures to illustrate convergence where the definition is via upper and lower sums (Sections D through I).
In Sections J and K, graphic procedures are provided to illustrate the Trapezoidal Rule and Simpson's Rule.
In Sections L, M and N, graphic procedures are provided to illustrate the use of the definite integral to represent area between curves, arc length, and area in polar coordinates.
A. Animating the
leftbox
,
rightbox
and
middlebox
commands
>
restart;
Animations are simple to create. On the whole, creating a sequence of plots, and then displaying them using the
display
command with the option
insequence=true
works better than using the animate command, which is rather restricted in its applications. For example suppose we want to illustrate convergence of Riemann sums approximating
. working from the graphics available in the
student
library, we have the choice of using left end-points, right end-points or midpoints over a uniform partition of the interval
. Begin by deciding on a sequence of values for
n
, the number of subintervals in the partition. There is little point in using a value of
n
greater than 100, as the rectangles will merge into a solid blob by then. For what follows, we will assume a standard sequence of values of
n
, namely:
>
L:=[5, 10, 20, 40, 60, 80, 100]:
Users may modify this to suit any particular needs. Using the
leftbox
command results in rectangular approximations where the function value at the left endpoint of each interval determines the height of the approximating rectangle over the subinterval. We use the integral noted above for illustration. Any continuous function on any finite interval may be used here.
>
f:=x->sin(x^2):a:=0:b:=sqrt(2*Pi):
We next construct appropriate plots for each
n
in
L.
>
with(student,leftbox):
for n in L
do Lb[n]:=leftbox(f(x),x=a..b,n):
od:
Finally, we use the
display
command with the option
insequence=true
to animate the sequential viewing of these plots.
>
with(plots,display):
display([seq(Lb[n],n=L)],insequence=true);
Some customization of the
leftbox
command is possible. Details are given on the relevant help page.
The related
rightbox
command uses right end-points, while the
middlebox
command uses midpoints of each sub-interval.
B. Riemann Sums with Uniform Partitions
>
restart;
Warning, the name changecoords has been redefined
Most standard textbooks define the definite integral as a limit of Riemann sums, so it is nice to have a procedure to visualize such sums. Here we consider the case where the partition is still uniform, but where sample points are chosen randomly in each subinterval. This is the definition used in, for instance
Stewart
, 4th Edition.
Consider a function
f
which is continuous for
x
in [
a
,
b
], with
a
<
b
. Partition the interval [
a
,
b
] into
n
sub-intervals, each of equal length
. In each sub-interval, select a sample point
, so that
. Then
In this section, we develop Maple tools to work with this definition both graphically and numerically, using a random process to select the sample points.
Rectangle Approximation Plot
Randbox
The inputs to the procedure are the function
f
, lower limit of integration
a
, upper limit
b
, and number of subintervals
n
in the uniform partition.
>
randbox:= proc(f,a,b,n)
local delta_x,i,j,k,B,C,A,x; #declares local variables
global x_star, randsum; #declares global variables
delta_x:=evalf((b-a)/n): #sets fineness of partition
with(plottools):
with(plots,display):
for i from 0 to n do
x[i]:=a+i*delta_x:od:i:='i': #names partition points
for i from 1 to n do
x_star[i]:=x[i-1]+rand()*10^(-12)*delta_x;od:i:='i': #selects sample points
randsum:=delta_x*sum(f(x_star[i]),i=1..n):i:='i': #sum of rectangle areas
for i from 1 to n do
B[i]:=POLYGONS([[x[i-1],0],[x[i-1],f(x_star[i])],[x[i],f(x_star[i])],
[x[i],0]]);od:i:='i': #builds rectangles
C:=PLOT(seq(B[i],i=1..n),COLOR(RGB,0,.8,.2)):
#plots rectangles
A:=plot(f(x),x=a..b,scaling=constrained,colour=red,thickness=2):
display({A,C},tickmarks=[3,3],labels=[``,``]); #combines plots
end:
A call to
randbox
now produces a visualization of one Riemann sum using a uniform partition with sample points randomly selected as noted above. A call to
randsum
after execution of randbox will produce the Riemann sum in decimal form. The command
[seq(x_star[i],i=1..n)];
after execution will display the sample points for this particular Riemann sum.
>
f:=x->sin(x^2):a:=0:b:=sqrt(2*Pi):
>
randbox(f,a,b,10);'randsum'=randsum;
`sample points`=[seq(x_star[j],j=1..10)];
Successive calls of the same
randbox
command will result in different plots, as the random numbers are newly generated at each call. Selection of the window display, instead of online, will allow direct comparison of such plots as separate windows for each are generated.
Technical Details
A procedure named
randbox
is defined with local variables, which revert to unassigned status after the procedure has been executed, and global variables, which remain assigned subsequently, and are thus retrievable at the user's option.
Following fairly standard notation, the width of one subinterval is called
, and the partition points
are generated by a do-loop. Sample points
are selected using the
rand
command, each call to which produces a random 12-digit integer, multiplied by
, to produce a random number in the interval (0, 1), and multiplied by
to produce a random positive number less than
. A do-loop calls this iteratively and adds the result to
to generate randomly and independently sample points
in each sub-interval. Since the variable
x_star
is declared as global, these sample points may be retrieved after executing the procedure if one wishes to see them.
The next do-loop used the
POLYGONS
command from the
plottools
library to build ingredients of the plot structure which will produce the
randbox
graphic. Each approximating rectangle is given its own label, and the
seq
command is used to put them all into an array, for which the
PLOT
command from the
plottools
library, not to be confused with the standard
plot
command, constructs a plot structure, in which a
COLOR
command, also from
plottools
, shades them green. This can be changed if desired; see the help page for
plot[structure]
for details.
The graph of the curve is plotted separately in a contrasting colour, and the two are assembled using
display
. If desired, tickmarks, axis labelling, etc, may be controlled by using the appropriate options in the
display
command.
Before ending the procedure, the signed sum of the rectangle areas is calculated and given the global variable name
randsum
. A call to
randsum
after the execution of
randbox
will then produce a floating-point value of the Riemann sum associated with the graphic produced.
Saving Procedures for future reference
Rather than reproducing the procedure each time, it may be saved using the command
>
save randbox, "c:\\randbox.m";
The procedure may be read in from such a file using the command
>
read "c:\\randbox.m";
Animating the
Randbox
Graphic
For a given function and interval, the convergence of these Riemann sums may be illustrated by animation as before. Note that each new call to the command will result in a different sequence of plots, as random numbers will be regenerated. Calling
sumlist
after execution of the animation will provide the sequence of signed rectangle sums from this particular execution.
>
L:=[5, 10, 20, 40, 60, 80, 100]:
>
sumlist:=[]:
for i in L do
RB[i]:=randbox(f,a,b,i):
sumlist:=[op(sumlist),randsum]:
od:
with(plots,display):
display([seq(RB[i],i=L)],insequence=true,scaling=constrained);
'sumlist'=sumlist;
Sums without Graphs
The
randsum
variable embedded in the
randbox
procedure can also be proceduralized without the graphics, so that numerical convergence can be observed beyond the limits of visual representation. To avoid conflict, when used on its own the procedure will be named
randSum
. Note that to get the Riemann sum associated with a particular graphic,
randsum
must be called after execution of the graphic. An independent execution of
randSum
will generate new sample points, and will thus not be the same Riemann sum as in the graphic.
The procedure is simply that as used before, omitting the graphics. To avoid accumulated roundoff error, precision is doubled. Maple automatically restores precision to that previous to execution when the procedure is exited.
>
restart;
>
randSum:= proc(f,a,b,n)
local delta_x,i,j,k,B,C,A,x; #declares local variables
global x_star; #declares global variable
Digits:=2*Digits:
delta_x:=evalf((b-a)/n): #sets fineness of partition
for i from 0 to n do
x[i]:=a+i*delta_x:od:i:='i': #names partition points
for i from 1 to n do
x_star[i]:=x[i-1]+rand()*10^(-12)*delta_x;od:i:='i': #selects sample points
delta_x*sum(f(x_star[i]),i=1..n); #sum of rectangle areas
end:
>
f:=x->sin(x^2):a:=0:b:=sqrt(2*Pi):
randSum(f,a,b,1000);Digits;
>
Int(f(x),x=a..b)=evalf(int(f(x),x=a..b));
C. Riemann Sums with Randomly Constructed Partitions
Some textbooks, such as
Stewart
, 3rd Edition, define the definite integral as a limit of Riemann sums over non-uniform partitions, as the length of the longest sub-interval, the norm of the partition, approaches 0. For a visualization of such situations, the following procedure constructs rectangle plots where the partition is randomly constructed, subject to a bound
on its norm, and sample points are randomly chosen in each sub-interval thus generated.
Rectangle Approximation Plot "Riembox"
The input to the procedure, caled
Riembox
, is the function, the limits of integration, and a positive number
to be taken as a bound on the norm. The procedure constructs a sequence of random positive numbers less than the norm, and uses them to determine a sequence of partition points from left to right. As soon as the most recently constructed partition point is past the right limit of integration, the partition is complete. There is no control over the number of sub-intervals, but the reality that sub-interval lengths are at least the norm divided by
ensures that the partition formation phase actually halts. The procedure then independently selects random sample points in each sub-interval, and constructs a rectangle plot as before. The global variable
Riemsum
reports the resulting sum in decimal form, while the global variable
normP
reports the actual norm, bounded by the number determined at the outset. The procedure also allows global variables to extract other information about the partition.
>
restart;
>
Riembox:= proc(f,a,b,delta) #names the procedure
local x,i,n,r,A,B,C; #declares local variables
global cellcount,P,delta_x,cell_sizes,normP,x_star,Riemsum; #declares global variable,
with(plottools):
x[0]:=evalf(a):P:=[x[0]]: #a is the first point in partition
for i from 0 while evalf(b)-x[i]>delta do
x[i+1]:=x[i]+delta*10^(-12)*rand(); #partition points chosen randomly
P:=[op(P),x[i+1]];
od:
cellcount:=nops(P); #number of sub-intervals
n:=cellcount:
x[n]:=evalf(b): #b is the last point
P:=[op(P),x[n]]: #The partition as a sequence
for i from 1 to n do
delta_x[i]:=x[i]-x[i-1]:od: #length of each sub-interval
cell_sizes:=[seq(delta_x[i],i=1..n)]: #list of interval lengths
normP:=max(seq(x[i]-x[i-1],i=1..n)); #norm of the partition
for i from 1 to n do
x_star[i]:=x[i-1]+10^(-12)*rand()*delta_x[i];
od: #Sample points chosen randomly in each sub-interval
Riemsum:=sum(f(x_star[k])*delta_x[k],k=1..n); #Riemann sum calculated
for i from 1 to n do
B[i]:=POLYGONS([[x[i-1],0],[x[i-1],f(x_star[i])],[x[i],f(x_star[i])],[x[i],0]]);od: #rectangle plot constructed
C:=PLOT(seq(B[i],i=1..n),COLOR(RGB,0,.8,.2)): #curve plotted
A:=plot(f(x),x=a..b,scaling=constrained,colour=red,thickness=2):
plots[display]({A,C}); #plots combined
end:
We apply this procedure to our (now standard) example. Note again that successive executions of this command produce new partitions and new sample points.
>
f:=x->sin(x^2):a:=0:b:=sqrt(2*Pi):
>
Riembox(f,a,b,0.3);`cell count`=cellcount;`cell sizes`=cell_sizes;
`norm of the partition`=normP;`Riemann sum`=Riemsum;
Animating the Riembox Plot
We can now animate this procedure. We need a list of partition norms approaching 0 through positive values this time. Some variables are delared as global, so numerical information about the partitions and sums can be extracted subsequently.
>
List:=[0.5,0.4,0.3,0.2,0.1,0.05]:
sumlist:=[]:normlist:=[]:
for delta in List
do Riem[delta]:=Riembox(f,a,b,delta):
sumlist:=[op(sumlist),Riemsum]:
normlist:=[op(normlist),normP]:
od:
plots[display]([seq(Riem[delta],delta=List)],insequence=true);
'normlist'=normlist;
'sumlist'=sumlist;
Sums Without Plots
To follow convergence numerically beyond what is visually observable, an independent
RiemSum
command is provided, which omits the plot structure. Other information about the partition can also be extracted as shown.
>
restart:
RiemSum:= proc(f,a,b,delta) #names the procedure
local x,i,k,n,r,A,B,C; #declares local variables
global cellcount,P,delta_x,cell_sizes,normP,x_star,sample_points; #declares global variable
Digits:=2*Digits:
x[0]:=evalf(a):P:=[x[0]]: #a is the first point in partition
for i from 0 while evalf(b)-x[i]>delta do
x[i+1]:=x[i]+delta*10^(-12)*rand(); #partition points chosen randomly
P:=[op(P),x[i+1]];
od:
cellcount:=nops(P); #number of sub-intervals
n:=cellcount:
x[n]:=evalf(b): #b is the last point
P:=[op(P),x[n]]: #The partition as a sequence
for i from 1 to n do
delta_x[i]:=x[i]-x[i-1]:od: #length of each sub-interval
cell_sizes:=[seq(delta_x[i],i=1..n)]: #list of interval lengths
normP:=max(seq(x[i]-x[i-1],i=1..n)); #norm of the partition
for i from 1 to n do
x_star[i]:=x[i-1]+10^(-12)*rand()*delta_x[i];
od: #Sample points chosen randomly in each sub-interval
sample_points:=[seq(x_star[i],i=1..n)];
sum(f(x_star[k])*delta_x[k],k=1..n); #Riemann sum calculated
end:
>
f:=x->sin(x^2):a:=0:b:=sqrt(2*Pi):
`Riemann Sum` = RiemSum(f,a,b,0.3);
`Number of Sub-intervals` = cellcount;
`Norm of the Partition` = normP;
`Partition Points` = [P];
`Sub-Interval Sizes` = cell_sizes;
`Sample Points` = sample_points;
D. Upper Sums with Uniform Partitions
The theory of definite integrals often requires the use of upper and lower sums as an alternate to Riemann sums for definition purposes. By their nature these require more care in Maple to illustrate. Any use of the procedures which follow should be limited to functions which are differentiable on the chosen interval, and have a finite number of zeros there.
First we look at upper sums for uniform partitions. These are also Riemann sums where the sample point is that which produces the maximum function value on the sub-interval. If this is not at an end point, then it is at a point where the derivative is 0, and there are finitely many of these if the restrictions just mentioned are observed. Maple can then quickly select the maximum point and construct the required rectangle.
Plotting Upper Sums: The Upbox Command
The following procedure,
upbox
, constructs the plot illustrating the upper sum. Inputs are the function, the limits of integration, and the number of subintervals in the uniform partition. Calling
upsum
reports the numerical value of the upper sum.
>
restart;
>
upbox:=proc(f,a,b,n)
local a1, b1, k, x, i, j, J, halt, critpt, critset, M, A, B, C;
global upsum;
with(plottools):
a1:=evalf(a):b1:=evalf(b): #puts limits in decimal form
k:=(b1-a1)/n; #sets fineness of partition
for i from 0 to n do
x[i]:=a1+i*k:od: #names partition points
critpt:=[];halt:=1;
while halt=1 do if type(fsolve(D(f)(x),x,avoid={seq(x=critpt[i],i=1..nops(critpt))},x=a1..b1),numeric)=true then critpt:=sort([op(critpt),fsolve(D(f)(x),x,avoid={seq(x=critpt[i],i=1..nops(critpt))},x=a1..b1)]);
else critpt; halt:=0;
fi;
od; #the list of zeros of the derivative on the interval is found
for i from 1 to n do
critset[i]:={x[i-1],x[i]};
for j in critpt do
if x[i-1]< j and j < x[i]
then critset[i]:={op(critset[i]),j};
else fi; #critical points in each subinterval are found
M[i]:=max(op(map(f,critset[i])));
od; #maxima are found in each subinterval
B[i]:=POLYGONS([[x[i-1],0],[x[i-1],M[i]],[x[i],M[i]],[x[i],0]]);
#rectangles are constructed
od:
upsum:=k*sum(M[J],J=1..n): #sum is calculated
C:=PLOT(seq(B[i],i=1..n),COLOR(RGB,0,.8,.2)):
A:=plot(f(x),x=a1..b1,scaling=constrained,colour=red,thickness=2):
plots[display]({A,C});#combines rectangle plot with plot of function
end:
>
f:=x->sin(x^2):a:=0:b:=sqrt(2*Pi):n:=10:
>
upbox(f,a,b,n);'upsum'=upsum;
Technical Details
For upper sums we need to find the maximum value of the function in each cell of the partition. These exist as soon as the function is continuous; if, in addition, the function is differentiable, the maxima occur either at end points or where the derivative is 0. We thus begin by listing all the zeros of the derivative over the complete interval. This is done by repetitive use of the
fsolve
command on the derivative. To build a list of zeros, we use the
avoid
option in
fsolve
at each iteration. Each time
fsolve
finds another zero, its output is of type
numeric
, while when there are no further zeros, the output of
fsolve
echoes the command, and thus is not of type
numeric.
The
type
command can thus be used to formulate a halting condition on this part of the process. The resulting list of zeros of the derivative is named
critpt.
For each sub-interval, a set of critical points
critset[i]
is formed, consisting of the end points of the sub-interval together with any zeros of the derivative lying in the subinterval. A comparison of function values yields the required maximum on the subinterval. Rectangles are then plotted, and sums calculated.
Animating Upper Sum Plots
To see the convergence of upper sums, this procedure can now be placed in an animation, and the sequence of upper sums thus formed reported.
>
L:=[5, 10, 20,40,60,80,100]:
>
sumlist:=[]:
for i in L do
UB[i]:=upbox(f,a,b,i):
sumlist:=[op(sumlist),upsum]:
od:
plots[display]([seq(UB[i],i=L)],insequence=true,scaling=constrained);
'sumlist'=sumlist;
Upper Sums Without Plots
This procedure, called
UpSum
, allows us to pursue upper sums numerically.
>
restart:
upSum:=proc(f,a,b,n)
local a1, b1, k, x, i, j, J, halt, critpt, critset, M;
Digits:=2*Digits:
a1:=evalf(a):b1:=evalf(b): #puts limits in decimal form
k:=(b1-a1)/n; #sets fineness of partition
for i from 0 to n
do x[i]:=a1+i*k:od: #names partition points
critpt:=[];halt:=1;
while halt=1 do
if type(fsolve(D(f)(x),x,avoid={seq(x=critpt[i],i=1..nops(critpt))},x=a1..b1),numeric)=true
then critpt:=sort([op(critpt),fsolve(D(f)(x),x,avoid={seq(x=critpt[i],i=1..nops(critpt))},x=a1..b1)]);
else critpt; halt:=0;
fi;
od; #the list of zeros of the derivative on the interval is found
for i from 1 to n do
critset[i]:={x[i-1],x[i]};
for j in critpt do
if x[i-1]< j and j < x[i]
then critset[i]:={op(critset[i]),j};
else fi; #critical points in each subinterval are found
M[i]:=max(op(map(f,critset[i])));
od;od; #maxima are found in each subinterval
k*sum(M[J],J=1..n); #sum is calculated
end:
>
f:=x->sin(x^2):a:=0:b:=sqrt(2*Pi):n:=10:
`Upper Sum` = upSum(f,a,b,n);
E. Lower Sums with Uniform Partitions
Analogously, here is a procedure, called
lowbox
, for illustrating lower sums with a uniform partition. The related numerical command is then
lowSum
.
Plotting Lower Sums
The method is similar to that for upper sums, using minima instead of maxima.
>
restart;
>
lowbox:=proc(f,a,b,n)
local a1, b1, k, x, i, j, J, halt, critpt, critset, m, A, B, C;#local variables
global lowsum; #declares global variables
with(plottools):
a1:=evalf(a):b1:=evalf(b): #puts limits in decimal form
k:=(b1-a1)/n; #sets fineness of partition
for i from 0 to n do
x[i]:=a1+i*k:od: #names partition points
critpt:=[];halt:=1;
while halt=1 do if type(fsolve(D(f)(x),x,avoid={seq(x=critpt[i],i=1..nops(critpt))},x=a1..b1),numeric)=true then critpt:=sort([op(critpt),fsolve(D(f)(x),x,avoid={seq(x=critpt[i],i=1..nops(critpt))},x=a1..b1)]);
else critpt; halt:=0;
fi;
od; #the list of zeros of the derivative on the interval is found
for i from 1 to n do
critset[i]:={x[i-1],x[i]};
for j in critpt do
if x[i-1]< j and j < x[i]
then critset[i]:={op(critset[i]),j};
else fi; #critical points in each subinterval are found
m[i]:=min(op(map(f,critset[i])));
od; #maxima are found in each subinterval
B[i]:=POLYGONS([[x[i-1],0],[x[i-1],m[i]],[x[i],m[i]],[x[i],0]]);
#rectangles are constructed
od:
lowsum:=k*sum(m[J],J=1..n): #sum is calculated
C:=PLOT(seq(B[i],i=1..n),COLOR(RGB,0,.8,.2)):
A:=plot(f(x),x=a1..b1,scaling=constrained,colour=red,thickness=2):
plots[display]({A,C});#combines rectangle plot with plot of function
end:
>
f:=x->sin(x^2):a:=0:b:=sqrt(2*Pi):n:=10:
>
lowbox(f,a,b,n);'lowsum'=lowsum;
Animating Lower Sum Plots
Again this procedure may be put into an animation.
>
L:=[5, 10, 20,40,60,80,100]:
>
sumlist:=[]:
for i in L do
LB[i]:=lowbox(f,a,b,i):
sumlist:=[op(sumlist),lowsum]:
od:
plots[display]([seq(LB[i],i=L)],insequence=true,scaling=constrained);
'sumlist'=sumlist;
Lower Sums Without Plots
The
lowSum
command allows lower sums to be calculated to greater accuracy than graphical displays permit.
>
restart;
>
lowSum:=proc(f,a,b,n)
local a1, b1, k, x, i, j, J, halt, critpt, critset, m;#local variables
Digits:=2*Digits:
a1:=evalf(a):b1:=evalf(b): #puts limits in decimal form
k:=(b1-a1)/n; #sets fineness of partition
for i from 0 to n do
x[i]:=a1+i*k:od: #names partition points
critpt:=[];halt:=1;
while halt=1 do if type(fsolve(D(f)(x),x,avoid={seq(x=critpt[i],i=1..nops(critpt))},x=a1..b1),numeric)=true then critpt:=sort([op(critpt),fsolve(D(f)(x),x,avoid={seq(x=critpt[i],i=1..nops(critpt))},x=a1..b1)]);
else critpt; halt:=0;
fi;
od; #the list of zeros of the derivative on the interval is found
for i from 1 to n do
critset[i]:={x[i-1],x[i]};
for j in critpt do
if x[i-1]< j and j < x[i]
then critset[i]:={op(critset[i]),j};
else fi; #critical points in each subinterval are found
m[i]:=min(op(map(f,critset[i])));
od;od; #maxima are found in each subinterval
k*sum(m[J],J=1..n):
end:
>
f:=x->sin(x^2):a:=0:b:=sqrt(2*Pi):n:=10:
`Lower Sum` = lowSum(f,a,b,n);
F. Difference of Upper and Lower Sums with Uniform Partitions
The proof that continuity implies the existence of the definite integral usually involves showing that the difference between upper and lower sums converges to 0 as the fineness of the partition decreases.
Plotting the Difference of Upper and Lower Sums
For the sake of illustration, the following procedure, called
UMLbox
(for Upper Minus Lower), is provided.
>
restart;
>
UMLbox:=proc(f,a,b,n)
local a1, b1, k, x, i, j, J, halt, critpt, critset,m, M, A, B, C;
global UMLsum;
with(plottools):
a1:=evalf(a):b1:=evalf(b): #puts limits in decimal form
k:=(b1-a1)/n; #sets fineness of partition
for i from 0 to n do
x[i]:=a1+i*k:od: #names partition points
critpt:=[];halt:=1;
while halt=1 do if type(fsolve(D(f)(x),x,avoid={seq(x=critpt[i],i=1..nops(critpt))},x=a1..b1),numeric)=true then critpt:=sort([op(critpt),fsolve(D(f)(x),x,avoid={seq(x=critpt[i],i=1..nops(critpt))},x=a1..b1)]);
else critpt; halt:=0;
fi;
od; #the list of zeros of the derivative on the interval is found
for i from 1 to n do
critset[i]:={x[i-1],x[i]};
for j in critpt do
if x[i-1]< j and j < x[i]
then critset[i]:={op(critset[i]),j};
else fi; #critical points in each subinterval are found
M[i]:=max(op(map(f,critset[i])));
m[i]:=min(op(map(f,critset[i])));
od; #maxima and minima are found in each subinterval
B[i]:=POLYGONS([[x[i-1],m[i]],[x[i-1],M[i]],[x[i],M[i]],[x[i],m[i]],[x[i-1],m[i]]]);
#rectangles are constructed
od:
UMLsum:=k*sum(M[J]-m[J],J=1..n): #sum is calculated
C:=PLOT(seq(B[i],i=1..n),COLOR(RGB,0,.8,.2)):
A:=plot(f(x),x=a1..b1,scaling=constrained,colour=red,thickness=2):
plots[display]({A,C});#combines rectangle plot with plot of function
end:
>
f:=x->sin(x^2):a:=0:b:=sqrt(2*Pi):n:=10:
>
UMLbox(f,a,b,n);'UMLsum'=UMLsum;
Animating the Difference Plot
Now we put this procedure into an animation.
>
L:=[5, 10, 20,40,60,80,100]:
sumlist:=[]:
for i in L do
UML[i]:=UMLbox(f,a,b,i):
sumlist:=[op(sumlist),UMLsum]:
od:
plots[display]([seq(UML[i],i=L)],insequence=true,scaling=constrained);
'sumlist'=sumlist;
Differences Without Plots
This command calculates the difference of upper and lower sums without plotting them. Double precision is used to avoid accumulated roundoff error in long sums.
>
restart;
>
UMLSum:=proc(f,a,b,n)
local a1, b1, k, x, i, j, J, halt, critpt, critset,m, M;
Digits:=2*Digits:
a1:=evalf(a):b1:=evalf(b): #puts limits in decimal form
k:=(b1-a1)/n; #sets fineness of partition
for i from 0 to n do
x[i]:=a1+i*k:od: #names partition points
critpt:=[];halt:=1;
while halt=1 do if type(fsolve(D(f)(x),x,avoid={seq(x=critpt[i],i=1..nops(critpt))},x=a1..b1),numeric)=true then critpt:=sort([op(critpt),fsolve(D(f)(x),x,avoid={seq(x=critpt[i],i=1..nops(critpt))},x=a1..b1)]);
else critpt; halt:=0;
fi;
od; #the list of zeros of the derivative on the interval is found
for i from 1 to n do
critset[i]:={x[i-1],x[i]};
for j in critpt do
if x[i-1]< j and j < x[i]
then critset[i]:={op(critset[i]),j};
else fi; #critical points in each subinterval are found
M[i]:=max(op(map(f,critset[i])));
m[i]:=min(op(map(f,critset[i])));
od;od; #maxima and minima are found in each subinterval
k*sum(M[J]-m[J],J=1..n):
end:
>
f:=x->sin(x^2):a:=0:b:=sqrt(2*Pi):n:=10:
`Upper Sum minus Lower Sum` = UMLSum(f,a,b,n);
G. Upper Sums with Random Partitions
Instead of uniform partitions, random partitions as in Section C can be used with upper sums.
Plotting Upper Sums with Random Partitions
Here we construct a command
uprandbox
for this purpose.
>
restart;
>
uprandbox:=proc(f,a,b,delta)
local x,i,j,J,r,critpt,critset,halt,A,B,C,M; #declares local variables
global n,P,delta_x,cell_sizes,normP,uprandsum; #declares global variables
with(plottools):
x[0]:=evalf(a):P:=[x[0]]: #a is the first point in partition
for i from 0 while evalf(b)-x[i]>delta do
x[i+1]:=x[i]+delta*10^(-12)*rand(); #partition points chosen randomly
P:=[op(P),x[i+1]];
od:
n:=nops(P); #number of sub-intervals
x[n]:=evalf(b): #b is the last point
P:=[op(P),x[n]]: #The partition as a sequence
for i from 1 to n do
delta_x[i]:=x[i]-x[i-1]:od: #length of each sub-interval
cell_sizes:=[seq(delta_x[i],i=1..n)]: #list of interval lengths
normP:=max(seq(x[i]-x[i-1],i=1..n)); #norm of the partition
critpt:=[];halt:=1;
while halt=1 do if type(fsolve(D(f)(x),x,avoid={seq(x=critpt[i],i=1..nops(critpt))},x=a..b),numeric)=true then critpt:=sort([op(critpt),fsolve(D(f)(x),x,avoid={seq(x=critpt[i],i=1..nops(critpt))},x=a..b)]);
else critpt; halt:=0;
fi;
od; #the list of zeros of the derivative on the interval is found
for i from 1 to n do
critset[i]:={x[i-1],x[i]};
for j in critpt do
if x[i-1]< j and j < x[i]
then critset[i]:={op(critset[i]),j};
else fi; #critical points in each subinterval are found
M[i]:=max(op(map(f,critset[i])));
od; #maxima are found in each subinterval
B[i]:=POLYGONS([[x[i-1],0],[x[i-1],M[i]],[x[i],M[i]],[x[i],0]]);
#rectangles are constructed
od:
uprandsum:=sum(M[J]*delta_x[J],J=1..n): #sum is calculated
C:=PLOT(seq(B[i],i=1..n),COLOR(RGB,0,.8,.2)):
A:=plot(f(x),x=a..b,scaling=constrained,colour=red,thickness=2):
plots[display]({A,C});#combines rectangle plot with plot of function
end:
Now let's try it.
>
f:=x->sin(x^2):a:=0:b:=sqrt(2*Pi):
>
uprandbox(f,a,b,0.3);'upsum'=uprandsum;'norm'=normP;
Animating the Upper Sum Plot for Random Partitions
>
List:=[0.3,0.2,0.1,0.05,0.03,0.02]:
sumlist:=[]:
normlist:=[]:
for delta in List do
URB[delta]:=uprandbox(f,a,b,delta):
sumlist:=[op(sumlist),uprandsum]:
normlist:=[op(normlist),normP]:
od:
plots[display]([seq(URB[delta],delta=List)],insequence=true,scaling=constrained);'sumlist'=sumlist;'normlist'=normlist;
Upper Sums without Plots
The command for upper sums with random partitions is
uprandSum
. Precision is doubled to avoid accumulated round-off error in long sums.
>
restart;
>
uprandSum:=proc(f,a,b,delta)
local x,i,j,J,r,critpt,critset,halt,M; #declares local variables
global n,P,delta_x,cell_sizes,normP; #declares global variables
Digits:=2*Digits;
x[0]:=evalf(a):P:=[x[0]]: #a is the first point in partition
for i from 0 while evalf(b)-x[i]>delta do
x[i+1]:=x[i]+delta*10^(-12)*rand(); #partition points chosen randomly
P:=[op(P),x[i+1]];
od:
n:=nops(P); #number of sub-intervals
x[n]:=evalf(b): #b is the last point
P:=[op(P),x[n]]: #The partition as a sequence
for i from 1 to n do
delta_x[i]:=x[i]-x[i-1]:od: #length of each sub-interval
cell_sizes:=[seq(delta_x[i],i=1..n)]: #list of interval lengths
normP:=max(seq(x[i]-x[i-1],i=1..n)); #norm of the partition
critpt:=[];halt:=1;
while halt=1 do if type(fsolve(D(f)(x),x,avoid={seq(x=critpt[i],i=1..nops(critpt))},x=a..b),numeric)=true then critpt:=sort([op(critpt),fsolve(D(f)(x),x,avoid={seq(x=critpt[i],i=1..nops(critpt))},x=a..b)]);
else critpt; halt:=0;
fi;
od; #the list of zeros of the derivative on the interval is found
for i from 1 to n do
critset[i]:={x[i-1],x[i]};
for j in critpt do
if x[i-1]< j and j < x[i]
then critset[i]:={op(critset[i]),j};
else fi; #critical points in each subinterval are found
M[i]:=max(op(map(f,critset[i])));
od;od; #maxima are found in each subinterval
sum(M[J]*delta_x[J],J=1..n);
end:
>
f:=x->sin(x^2):a:=0:b:=sqrt(2*Pi):
`Upper Sum` = uprandSum(f,a,b,0.3);
H. Lower Sums with Random Partitions
Plotting Lower Sums with Random Partitions
>
restart;
>
lowrandbox:=proc(f,a,b,delta)
local x,i,j,J,r,critpt,critset,halt,A,B,C,m; #declares local variables
global n,P,delta_x,cell_sizes,normP,lowrandsum; #declares global variables
with(plottools):
x[0]:=evalf(a):P:=[x[0]]: #a is the first point in partition
for i from 0 while evalf(b)-x[i]>delta do
x[i+1]:=x[i]+delta*10^(-12)*rand(); #partition points chosen randomly
P:=[op(P),x[i+1]];
od:
n:=nops(P); #number of sub-intervals
x[n]:=evalf(b): #b is the last point
P:=[op(P),x[n]]: #The partition as a sequence
for i from 1 to n do
delta_x[i]:=x[i]-x[i-1]:od: #length of each sub-interval
cell_sizes:=[seq(delta_x[i],i=1..n)]: #list of interval lengths
normP:=max(seq(x[i]-x[i-1],i=1..n)); #norm of the partition
critpt:=[];halt:=1;
while halt=1 do if type(fsolve(D(f)(x),x,avoid={seq(x=critpt[i],i=1..nops(critpt))},x=a..b),numeric)=true then critpt:=sort([op(critpt),fsolve(D(f)(x),x,avoid={seq(x=critpt[i],i=1..nops(critpt))},x=a..b)]);
else critpt; halt:=0;
fi;
od; #the list of zeros of the derivative on the interval is found
for i from 1 to n do
critset[i]:={x[i-1],x[i]};
for j in critpt do
if x[i-1]< j and j < x[i]
then critset[i]:={op(critset[i]),j};
else fi; #critical points in each subinterval are found
m[i]:=min(op(map(f,critset[i])));
od; # minima are found in each subinterval
B[i]:=POLYGONS([[x[i-1],0],[x[i-1],m[i]],[x[i],m[i]],[x[i],0]]);
#rectangles are constructed
od:
lowrandsum:=sum(m[J]*delta_x[J],J=1..n): #sum is calculated
C:=PLOT(seq(B[i],i=1..n),COLOR(RGB,0,.8,.2)):
A:=plot(f(x),x=a..b,scaling=constrained,colour=red,thickness=2):
plots[display]({A,C});#combines rectangle plot with plot of function
end:
Applying this,
>
f:=x->sin(x^2):a:=0:b:=sqrt(2*Pi):
>
lowrandbox(f,a,b,0.3);'lowsum'=lowrandsum;'norm'=normP;
Animating the Plot
>
List:=[0.3,0.2,0.1,0.05,0.03,0.02]:
sumlist:=[]:
normlist:=[]:
for delta in List do
LRB[delta]:=lowrandbox(f,a,b,delta):
sumlist:=[op(sumlist),lowrandsum]:
normlist:=[op(normlist),normP]:
od:
plots[display]([seq(LRB[delta],delta=List)],insequence=true,scaling=constrained);'sumlist'=sumlist;'normlist'=normlist;
Lower Sums Without Plots
>
restart;
>
lowrandSum:=proc(f,a,b,delta)
local x,i,j,J,r,critpt,critset,halt,m; #declares local variables
global n,P,delta_x,cell_sizes,normP; #declares global variables
Digits:=2*Digits;
x[0]:=evalf(a):P:=[x[0]]: #a is the first point in partition
for i from 0 while evalf(b)-x[i]>delta do
x[i+1]:=x[i]+delta*10^(-12)*rand(); #partition points chosen randomly
P:=[op(P),x[i+1]];
od:
n:=nops(P); #number of sub-intervals
x[n]:=evalf(b): #b is the last point
P:=[op(P),x[n]]: #The partition as a sequence
for i from 1 to n do
delta_x[i]:=x[i]-x[i-1]:od: #length of each sub-interval
cell_sizes:=[seq(delta_x[i],i=1..n)]: #list of interval lengths
normP:=max(seq(x[i]-x[i-1],i=1..n)); #norm of the partition
critpt:=[];halt:=1;
while halt=1 do if type(fsolve(D(f)(x),x,avoid={seq(x=critpt[i],i=1..nops(critpt))},x=a..b),numeric)=true then critpt:=sort([op(critpt),fsolve(D(f)(x),x,avoid={seq(x=critpt[i],i=1..nops(critpt))},x=a..b)]);
else critpt; halt:=0;
fi;
od; #the list of zeros of the derivative on the interval is found
for i from 1 to n do
critset[i]:={x[i-1],x[i]};
for j in critpt do
if x[i-1]< j and j < x[i]
then critset[i]:={op(critset[i]),j};
else fi; #critical points in each subinterval are found
m[i]:=min(op(map(f,critset[i])));
od;od; # minima are found in each subinterval
sum(m[J]*delta_x[J],J=1..n);
end:
>
f:=x->sin(x^2):a:=0:b:=sqrt(2*Pi):delta:=0.3:
`Lower Sum` = lowrandSum(f,a,b,delta);
I. Difference of Upper and Lower Sums with Random Partitions
Plotting the Difference of Upper and Lower Sums with Random Partitions
>
restart;
>
UMLrandbox:=proc(f,a,b,delta)
local x,i,j,J,r,critpt,critset,halt,A,B,C,m,M; #declares local variables
global n,P,delta_x,cell_sizes,normP,UMLrandsum; #declares global variables
with(plottools):
x[0]:=evalf(a):P:=[x[0]]: #a is the first point in partition
for i from 0 while evalf(b)-x[i]>delta do
x[i+1]:=x[i]+delta*10^(-12)*rand(); #partition points chosen randomly
P:=[op(P),x[i+1]];
od:
n:=nops(P); #number of sub-intervals
x[n]:=evalf(b): #b is the last point
P:=[op(P),x[n]]: #The partition as a sequence
for i from 1 to n do
delta_x[i]:=x[i]-x[i-1]:od: #length of each sub-interval
cell_sizes:=[seq(delta_x[i],i=1..n)]: #list of interval lengths
normP:=max(seq(x[i]-x[i-1],i=1..n)); #norm of the partition
critpt:=[];halt:=1;
while halt=1 do if type(fsolve(D(f)(x),x,avoid={seq(x=critpt[i],i=1..nops(critpt))},x=a..b),numeric)=true then critpt:=sort([op(critpt),fsolve(D(f)(x),x,avoid={seq(x=critpt[i],i=1..nops(critpt))},x=a..b)]);
else critpt; halt:=0;
fi;
od; #the list of zeros of the derivative on the interval is found
for i from 1 to n do
critset[i]:={x[i-1],x[i]};
for j in critpt do
if x[i-1]< j and j < x[i]
then critset[i]:={op(critset[i]),j};
else fi; #critical points in each subinterval are found
m[i]:=min(op(map(f,critset[i])));
M[i]:=max(op(map(f,critset[i])));
od; # maxima and minima are found in each subinterval
B[i]:=POLYGONS([[x[i-1],m[i]],[x[i],m[i]],[x[i],M[i]],[x[i-1],M[i]]]);
#rectangles are constructed
od:
UMLrandsum:=sum((M[J]-m[J])*delta_x[J],J=1..n): #sum is calculated
C:=PLOT(seq(B[i],i=1..n),COLOR(RGB,0,.8,.2)):
A:=plot(f(x),x=a..b,scaling=constrained,colour=red,thickness=2):
plots[display]({A,C});#combines rectangle plot with plot of function
end:
>
f:=x->sin(x^2):a:=0:b:=sqrt(2*Pi):
>
UMLrandbox(f,a,b,0.3);'UMLsum'=UMLrandsum;'norm'=normP;`number of cells`=n;
Animating the Plot
>
List:=[0.3,0.2,0.1,0.05,0.03,0.02]:
sumlist:=[]:
normlist:=[]:
for delta in List do
UML[delta]:=UMLrandbox(f,a,b,delta):
sumlist:=[op(sumlist),UMLrandsum]:
normlist:=[op(normlist),normP]:
od:
plots[display]([seq(UML[delta],delta=List)],insequence=true,scaling=constrained);'sumlist'=sumlist;'normlist'=normlist;
Difference of Sums Without Plotting
>
restart;
>
UMLrandSum:=proc(f,a,b,delta)
local x,i,j,J,r,critpt,critset,halt,m,M; #declares local variables
global n,P,delta_x,cell_sizes,normP; #declares global variables
Digits:=2*Digits;
x[0]:=evalf(a):P:=[x[0]]: #a is the first point in partition
for i from 0 while evalf(b)-x[i]>delta do
x[i+1]:=x[i]+delta*10^(-12)*rand(); #partition points chosen randomly
P:=[op(P),x[i+1]];
od:
n:=nops(P); #number of sub-intervals
x[n]:=evalf(b): #b is the last point
P:=[op(P),x[n]]: #The partition as a sequence
for i from 1 to n do
delta_x[i]:=x[i]-x[i-1]:od: #length of each sub-interval
cell_sizes:=[seq(delta_x[i],i=1..n)]: #list of interval lengths
normP:=max(seq(x[i]-x[i-1],i=1..n)); #norm of the partition
critpt:=[];halt:=1;
while halt=1 do if type(fsolve(D(f)(x),x,avoid={seq(x=critpt[i],i=1..nops(critpt))},x=a..b),numeric)=true then critpt:=sort([op(critpt),fsolve(D(f)(x),x,avoid={seq(x=critpt[i],i=1..nops(critpt))},x=a..b)]);
else critpt; halt:=0;
fi;
od; #the list of zeros of the derivative on the interval is found
for i from 1 to n do
critset[i]:={x[i-1],x[i]};
for j in critpt do
if x[i-1]< j and j < x[i]
then critset[i]:={op(critset[i]),j};
else fi; #critical points in each subinterval are found
m[i]:=min(op(map(f,critset[i])));
M[i]:=max(op(map(f,critset[i])));
od;od; # maxima and minima are found in each subinterval
sum((M[J]-m[J])*delta_x[J],J=1..n);
end:
>
f:=x->sin(x^2):a:=0:b:=sqrt(2*Pi):delta:=0.3:
`Difference of Sums` = UMLrandSum(f,a,b,delta);
J. Trapezoidal Rule
Maple's Student library includes a
trapezoid
command to calculate the approximation to a definite integral using the Trapezoidal Rule. Here is a plot which illustrates the situation.
Plotting the Trapezoidal Rule
For best results, the function chosen here should be positive throughout the interval selected. The command
trapbox
constructs a plot showing the trapezoidal approximation.
>
restart;
>
trapbox:=proc(f,a,b,n)
local c,h,i,pol,pl; #local variables
with(plottools):
h:=(b-a)/n:
for i from 1 to n do pol[i]:=polygon([[a+(i-1)*h,0],[a+(i-1)*h,f(a+(i-1)*h)],[a+i*h,f(a+i*h)],[a+i*h,0]],colour=green):od:
pl[n]:=plots[display](PLOT(seq(pol[i],i=1..n)),scaling=constrained,tickmarks=[3,3]):
c:=plot(f(x),x=a..b,colour=red,thickness=2):
plots[display]({c,pl[n]},scaling=constrained);
end:
>
f:=x->1+sin(x^2):a:=0:b:=sqrt(2*Pi):n:=8:
trapbox(f,a,b,n);
Trap[8] = evalf(student[trapezoid](f(x),x=a..b,n));
Technical Details
The plot of each trapezoid is constructed using the
polygon
command in the
plottools
library. These are then assembled and displayed along with the plot of the original curve.
Animating the Trapbox Plot
To see that increasing the number of subdivisions improves the accuracy of the trapezoidal approximation we can animate the
trapbox
plots.
>
L:=[5, 10, 20, 40, 60]:
>
for n in L do
A[n]:=trapbox(f,a,b,n):
od:
plots[display]([seq(A[n],n=L)],insequence=true);
K. Simpson's Rule
Maple's Student library contains the
simpson
command to calculate the Simpson's Rule approximation to a definite integral. Here we construct a command
simpsonplot
to show the area below the approximating parabolas. Note that
n
subdivisions produce
approximating parabolas.
Plotting Simpson's Rule Approximations
For this procedure, it is best to select a function which is positive over the selected interval. The number of sub-intervals selected must be even.
>
restart;
>
simpsonplot:=proc(f,a,b,n)
local h,k,A,B,C,g,PAR,POL,P,P1,V,Vert,W,WHITE; #local variables
with(plottools):
h:=evalf((b-a)/n): #length of one subdivision
WHITE:=[]:
for k from 0 by 2 to n-2 do
A[k]:=1/2*(-2*f(a+k*h+h)+f(a+k*h)+f(a+k*h+2*h))/(h^2):
B[k]:=-1/2*(2*a*f(a+k*h+2*h)-4*a*f(a+k*h+h)+2*f(a+k*h)*a+3*f(a+k*h)*h-4*f(a+k*h+h)*h+2*f(a+k*h+2*h)*k*h+2*f(a+k*h)*k*h-4*f(a+k*h+h)*k*h+f(a+k*h+2*h)*h)/(h^2):
C[k]:=1/2*(-2*f(a+k*h+h)*a^2+f(a+k*h)*a^2+f(a+k*h+2*h)*a^2+3*f(a+k*h)*h*a-4*a*f(a+k*h+h)*h+2*f(a+k*h+2*h)*a*k*h+2*f(a+k*h)*a*k*h-4*f(a+k*h+h)*a*k*h+a*f(a+k*h+2*h)*h+f(a+k*h)*k^2*h^2+f(a+k*h+2*h)*k^2*h^2+f(a+k*h+2*h)*k*h^2+3*f(a+k*h)*k*h^2+2*f(a+k*h)*h^2-2*f(a+k*h+h)*k^2*h^2-4*f(a+k*h+h)*k*h^2)/(h^2): #coefficients of parabolic approximation
g[k]:=x->A[k]*x^2+B[k]*x+C[k]; #equation of parabolic approximation
PAR[k]:=plot(g[k](x),x=a+k*h..a+(k+2)*h):
POL[k]:=PLOT(polygon([op(op([1,1],PAR[k])),[a+(k+2)*h,0],[a+k*h,0]],colour=green));
V[k]:=plot([[a+(k+1)*h,0],[a+(k+1)*h,f(a+(k+1)*h)]],colour=black):
if A[k]>0 then
W[k]:=PLOT(polygon([op(op([1,1],PAR[k]))],colour=white,style=patchnogrid));
WHITE:=[op(WHITE),W[k]];
else fi;
od:
P:=plots[display]([seq(POL[2*k],k=0..(n-2)/2)]):
Vert:=plots[display]([seq(V[2*k],k=0..(n-2)/2)]):
P1:=plot(f(x),x=a..b,thickness=3,colour=red):
plots[display]([op(WHITE),P,P1,Vert],scaling=constrained);
end:
>
f:=x->1+sin(x^2):a:=0:b:=sqrt(2*Pi):n:=6:
simpsonplot(f,a,b,n);
S[4]=evalf(student[simpson](f(x),x=a..b,n));
Technical Details
The challenge here is to plot not just the approximating parabolas, but the area below them. After setting up the partition, parabolas are fitted to each triple of points on the curve. To do this, one takes the points
,
and
, where
h
is the sub-interval width and
k
the index of the sub-interval pair, and fits a parabola to them. The parabola then has equation
. The equations for
,
and
in terms of
f
,
a
,
h
and
k
were worked out using Maple separately.
A plot structure for each parabola over the appropriate sub-interval pair is next constructed; these are designated
in the procedure. The sequence of points used for the parabola plot is then extracted from each plot structure by careful use of the
op
command, and joined to the end points of the sub-interval pair on the
x
axis, inside a
polygon
command from
plottools.
What Maple actually plots is the convex hull of this polygon, so for those parabolic sec\ments which are concave upwards, we superimpose another plot in white
. Using the
PLOT
command from
plottools
, each region below a parabolic approximation is labelled
in the procedure. The sequence of these is then displayed along with the curve of the function and additional vertical lines to display the partition..
Animating Plots for Simpson's Rule
An increasing sequence of even values of
n
must be used here.
>
L:=[4,8,16,32]:
f:=x->1+sin(x^2):a:=0:b:=sqrt(2*Pi):
>
for n in L do
A[n]:=simpsonplot(f,a,b,n):
od:
plots[display]([seq(A[n],n=L)],insequence=true,scaling=constrained);
L. Area Between Curves
Here we present a command,
betweenbox
, for plotting rectangular approximations to areas between curves, using midpoints of regular partitions. Calculation of the corresponding sums can then be done by applying the
middlesum
command from the Student library.
Plotting Rectangular Approximations to Areas Between Curves
>
restart;
>
betweenbox:=proc(f,g,a,b,n)
local a1,b1,k,i,j,P,B,A;
with(plottools):
a1:=evalf(a):b1:=evalf(b):
k:=(b1-a1)/n:
for i in [seq(j,j=0..n-1)] do
P[i,n]:=POLYGONS([[a1+k*i,g(a1+k/2+k*i)],[a1+k*i,f(a1+k/2+k*i)],[a1+k*(i+1),f(a1+k/2+k*i)],[a1+k*(i+1),g(a1+k/2+k*i)]]);
od:
B[n]:=PLOT(seq(P[i,n],i=0..n-1),COLOR(RGB,0,.8,.2)):
A:=plot({f(x),g(x),[[a1,f(a1)],[a1,g(a1)]],[[b1,f(b1)],[b1,g(b1)]]},x=a1..b1,colour=red,scaling=constrained,thickness=2,labels=[" "," "]):
plots[display]({A,B[n]});end:
>
f:=x->sin(x):g:=x->cos(x):a:=Pi/4:b:=5*Pi/4:n:=10:
betweenbox(f,g,a,b,n);
`Rectangle Area`=evalf(student[middlesum](f(x)-g(x),x=a..b,n));
Technical Details
Using the
polygons
command in the
plottools
library, individual rectangle plot structures are constructed and assembled. These are then plotted using
display
.
Animating the Plots
>
L:=[5, 10, 20, 40, 60, 80, 100]:
>
for n in L
do BB[n]:=betweenbox(f,g,a,b,n):
od:
plots[display]([seq(BB[n],n=L)],insequence=true);
M. Arc Length
Here we present a procedure to illustrate the approximation of arc length by a sequence of straight line segments.
Plotting Arc Length Approximations
This procedure plots the curve together with approximating line segments, using a uniform partition. The length of the approximating segments can be extracted using the global variable
arcsum
.
>
restart;
>
arcplot:=proc(f,a,b,n)
local A,B,i,h;
global arcsum;
h:=(b-a)/n:
A:=plot(f(x),x=a..b,colour=red,thickness=2):
B:=plot([seq([a+i*h,f(a+i*h)],i=0..n)],colour=black):
arcsum:=evalf(sum(sqrt(h^2+(f(a+(i+1)*h)-f(a+i*h))^2),i=0..n-1)):
plots[display]({A,B});
end:
>
f:=x->sin(x^2):a:=0:b:=sqrt(2*Pi):n:=5:
arcplot(f,a,b,n);
`Line Approximation` = arcsum;
Animating the Convergence of Approximations
The convergence of the approximations can be illustrated using an increasing sequence of values of
n
. The corresponding numerical approximations can be generated at the same time.
>
L:=[5,10,15,20]:
for j in L do
C[j]:=arcplot(f,a,b,j):
S[j]:=arcsum:
od:
sumlist:=[seq(S[j],j=L)]:
plots[display]([seq(C[j],j=L)],insequence=true);
'Sumlist'=sumlist;
Approximating Sums Without Plots
If it is desired to follow the convergence of the approximating sums beyond the visual range, this command will allow them to be generated. Precision is doubled to avoid accumulated roundoff error in longer sums
>
restart;
>
arcSum:=proc(f,a,b,n)
local i,h;
Digits:=2*Digits;
h:=(b-a)/n:
evalf(sum(sqrt(h^2+(f(a+(i+1)*h)-f(a+i*h))^2),i=0..n-1)):
end:
>
f:=x->sin(x^2):a:=0:b:=sqrt(2*Pi):n:=5:
arcSum(f,a,b,n);
N. Area in Polar Coordinates
When setting up a definite integral to represent area in polar coordinates, the interval is a central angle, and the approximating rectangles become sectors of circles. The graphic command
polarbox
illustrates the procedure.
Plotting Approximations to Area in Polar Coordinates
Here we assume that
r
is given as a function of
, and is positive for the
interval [
a
,
b
]. The partition of this interval is uniform, and sample points are midpoints.
>
restart;
>
polarbox:=proc (r, a, b, n)
local a1, b1, k, i,j, P, A, B;
global sectorsum;
with(plottools,pieslice);
a1 := evalf(a); b1 := evalf(b); k := (b1-a1)/n;
for i to n do
P[i] := pieslice([0, 0],r(a1-1/2*k+k*i),a1+k*(i-1) .. a1+k*i,colour = COLOUR(RGB,0,.7,.3));
B[i]:=evalf(r(a1-1/2*k+k*i)^2*k/2); od:
A := plot([r(theta), theta, theta = a .. b],coords = polar,scaling = constrained,colour = red,thickness = 2):
sectorsum:=sum(B[j],j=1..n):
plots[display]({seq(P[i],i = 1 .. n), A},scaling = constrained);
end:
>
r:=theta->4*(1-cos(theta)):a:=0:b:=Pi:n:=6:
polarbox(r,0,Pi,6);
`Sector Sum` = sectorsum;
Technical Details
After setting up the partition, the
pieslice
command in the
plottools
library allows the plot structure for each approximating sector to be constructed. These are then assembled and plotted using
display
.
Animating the Plot
>
L:=[5, 10, 20, 40, 60, 80, 100]:
>
for n in L do
PP[n]:=polarbox(r,a,b,n):
od:
plots[display]([seq(PP[n],n=L)],insequence=true);
Approximating Sums Without Plots
To follow the convergence of an approximating sum beyond the visual range, the command
sectorSum
is provided. To avoid accumulated round-off error, precision is doubled.
>
restart;
>
sectorSum:=proc (r, a, b, n)
local a1, b1, k, i,j, P, B;
Digits:=2*Digits;
a1 := evalf(a); b1 := evalf(b); k := (b1-a1)/n;
for i to n do
B[i]:=evalf(r(a1-1/2*k+k*i)^2*k/2); od:
sum(B[j],j=1..n):
end:
>
r:=theta->4*(1-cos(theta)):a:=0:b:=Pi:n:=6:
sectorSum(r,a,b,n);