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
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.
> 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.
> 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.
> 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.
> 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.
> 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.
> 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;
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.
> 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:
> lowbox(f,a,b,n);'lowsum'=lowsum;
Animating Lower Sum Plots
Again this procedure may be put into an animation.
> 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.
> 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.
> 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:
> 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.
> 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.
> 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.
> 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.
> 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
> 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,
> 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;
> 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
> 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:
> UMLrandbox(f,a,b,0.3);'UMLsum'=UMLrandsum;'norm'=normP;`number of cells`=n;
> 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
> 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.
> 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));
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.
> 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));
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
> 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));
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
> 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 .
> 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
> 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.
> 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;
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 .
> for n in L do PP[n]:=polarbox(r,a,b,n): od: plots[display]([seq(PP[n],n=L)],insequence=true);
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.
> 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);
