Application Center - Maplesoft

Animations for the calculus classroom: definite integral approximations

You can switch back to the summary page by clicking here.

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

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);