Calculus Projects with Maple

30. Iterated Integrals

piecewise	plot3d	Doubleint
Int	int	value
rand	seq	union
sort	convert	?

In addition, you may find it useful to be able to use the list option of convert.

Jump to section:

1. Background
2. Extension
3. Background The double integral of a function f(x,y) over finite region R in the plane can be considered a geometric or an analytic object. As a geometric object, it is the (signed) volume of region in three space: that cylinder in three space bounded by R and the projection of R onto the surface z=f(x,y) . As an analytic object, the double integral is the limiting value of a collection of Riemann sums. Each of these Riemann sums is formed by tiling (nearly) R with rectangles, then adding the products of the area of each tile with the value of f at some point (x,y) chosen interior to that tile. These conceptualizations are intuitively appealing (as long as we don't have to evaluate the Riemann sums), complementary (one gives a nice picture, but is no help at computation, the other is purely computational), and clearly related (volume is area times height). We evaluate these geometric/analytic objects, double integrals, using the purely algebraic device of a twice iterated integral. Typically, the technical difficulties of finding an algebraic description of R , of choosing an order for the integration, "first x , then y ," or "first y , then x," etc., distract us from the fact that there is no obvious assurance that our algebraic technique reflects the geometric/analytic reality. We have, from the time we first studied algebra, become so used to seeing geometry reflected in our algebra that we seldom question the efficacy of algebra to describe geometry. In this project we will examine a fairly simple function for which the algebraic technique of twice iterated integrals really does fail to capture the geometric reality of the double integral. Such functions are not hard to construct, but they are hard to recognize. We will also examine a very similar function for which the iterated integrals do capture the double integral. Lest you worry that iterated integration is totally useless, we hasten to point out that the double integral of any smooth function over any region with a nice boundary is faithfully captured by either, and hence both, of the iterated integrals. This result is known as Fubini's Theorem (Guido Fubini (1879--1943) was an Italian mathematician. His theorem, published in 1907, is more general, not requiring f to be continuous.).

   Project Report

   We'll consider the function f defined on the unit square R by

   

   This bad function (some describe it as pathological, but that seems harsh to us) is of necessity discontinuous, so we describe it as a piecewise function. Maple has a command that facilitates the construction of piecewise functions, it is called---not surprisingly---piecewise. Define f(x,y) by: (Release 3 users, jump to Alternate Function Definition.)

   f := (x,y) -> piecewise(
   	0<x and x<y and y<1, 1/y^2,
   	0<y and y<x and x<1, -1/x^2);
   

   Actually, f is defined to be zero everywhere else---particularly on the diagonal x=y and on the edges of the unit square---as piecewise assumes. However, we get a better graph if we don't let Maple use the default 0. Specify an unevaluated result by using `F(x,y)` as the "otherwise" clause.

   F := (x,y) -> piecewise(
   	0<x and x<y and y<1, 1/y^2,
   	0<y and y<x and x<1, -1/x^2,
   	`F(x,y)`):

   Now to obtain the image, enter:

   plot3d(F, 0..1, 0..1, view=-10..10, labels=[X,Y,Z], axes=frame, numpoints=900);

   Alternate Function Definition

   There is a flaw in Release 3's piecewise; compound conditions aren't accepted. It will be easier to use a compound if-then. Replace f's definition above with:

   f := (x,y) -> if 0<x and x<y and y<1 then 1/y^2
   	elif 0<y and y<x and x<1 then -1/x^2
   	else 0 fi:

   Also replace F's definition with:

   F := (x,y) -> if 0<x and x<y and y<1 then 1/y^2
   	elif 0<y and y<x and x<1 then -1/x^2
   	else `F(x,y)` fi:

   There are several obvious items to note. The graph shows two separate sheets; f is discontinuous as promised. The surface exits to + infinity and to - infinity as x and y approach the origin. Fubini's Theorem requires that the discontinuity be real---not just a gratuitous shift of a smooth surface. The graph has odd symmetry with respect to the diagonal line y=x . That is, the positive half of the graph is a reflection of the the negative half through the line y=x . The double integral ought evaluate to zero---if it exists, not all improper integrals converge. To see what the algebra does, look at both iterated integrals:
   • Integrate x first then y, denoted by ,
   • integrate y first then x, denoted by .
   There is a twice iterated integral command in the student package listed under the customary name of Doubleint.

   with(student):
   
   xFirst := Doubleint(f(x,y), x=0..1, y=0..1);
   yFirst := Doubleint(f(x,y), y=0..1, x=0..1);
   

   

   How you proceed from here depends on the size and speed of the machine you are using and on the version of Maple available. If you are running Maple V Release 4 and have a sleek, powerful mainframe with lots of memory (or if you merely have the time), ask for a numeric value of each of the integrals with evalf. If your resources are like those of the authors, more modest---we are using a Power Macintosh 7100/66---follow the logic through below. To evaluate xFirst:

   Int(1/y^2, x=0..y) + Int(-1/x^2, x=y..1);
   value(");

   

   1

   Int(", y=0..1);
   value(");

   

   1

   Interesting---not the value we were expecting.

   Question 1 : Explain why the above is a valid computation for xFirst and the following is a valid computation for yFirst.

   Now for yFirst.

   Int(-1/x^2, y=0..x) + Int(1/y^2, y=x..1);
   value(");

   

   Int(", x=0..1);
   value(");

   

   This is positively distressing! Not only do the separate computations of xFirst and yFirst fail to give the geometrically plausible value of 0, they don't give the same (wrong) answer. This conundrum must be due to the improper integrals. However, that is not the entire story; consider the function g that follows. (Release 3 users, jump to Alternate Definition.)

   g := (x,y) -> piecewise(
   	0<x and x<y and y<1, 1/y^(1/2),
   	0<y and y<x and x<1, -1/x^(1/2)):

   As before, to obtain a better graph use:

   G := (x,y) -> piecewise(
   	0<x and x<y and y<1, 1/y^(1/2),
   	0<y and y<x and x<1, -1/x^(1/2),
   	`G(x,y)`):
   
   plot3d(G, 0..1, 0..1, view=-10..10, labels = [X,Y,Z], axes=frame, numpoints=900);

   Alternate Definition

   Use compound if-then definitions for g and G paralleling those you used for f and F previously.

   The graph of g looks similar to f's with a less dramatic exit at the origin, but still going for infinite values. The twice iterated integrals parallel those of f.

   xFirstG := Int(1/y^(1/2), x=0..y) + Int(-1/x^(1/2), x=y..1);
   value(");

   

   Int(", y=0..1);
   value(");

   

   and

   yFirstG := Int(-1/x^(1/2), y=0..x) + Int(1/y^(1/2), y=x..1);
   value(");

   

   Int(", x=0..1);
   value(");

   

   0

   The iterated integrals each agree with the geometrically plausible 0. To gain insight into the difficulties with f , we return to Riemann sums and use a two variable variation of the Monte Carlo integration we've looked at in previous projects. This method---refined and rationalized---is actually one of the preferred methods for evaluating multiple integrals in more than two or three variables. We will do a small sampling here and ask you and your class to collect a larger data set for further analysis.

   _seed := 123456789:
   r := () -> rand()/10.0^12:
   This defines r as a random number generator that returns decimal numbers
   between 0 and 1.
   
   X := {seq(r(), k=1..99)} union {0, 1}:
   X := sort(convert(X, list));
   
   	
   
   
   Now we have a list of values to use for the vertical grid lines of our
   tiling and for the x -coordinates of the sample points within the tiles. We
   intend to use the odd items in this list for the grid lines and the
   even items for the x -coordinates. Do the same to obtain horizontal grid
   lines and y -coordinates. 
   Y := {seq(r(), k=0..100)} union {0, 1}:
   Y := sort(convert(Y, list)):
   
   Calculate the general term---one for each tile in our grid---of the Riemann sum.
   This term is the product of the function evaluated at the sample point in the tile
   with the area of the tile itself.
   
   generalTerm := f(X[2*n],Y[2*n])
   	* (X[2*n+1]-X[2*n-1]) * (Y[2*n+1]-Y[2*n-1]):
   s := sum(generalTerm, n=1..50);
   
   	
   s := -1.124350811
   
   
   Lets do this again.
   X := {seq(r(), k=1..99)} union {0, 1}:
   X := sort(convert(X,list)):
   Y := {seq(r(), k=0..100)} union {0, 1}:
   Y := sort(convert(Y, list)):
   generalTerm := f(X[2*n],Y[2*n])
   	* (X[2*n+1]-X[2*n-1]) * (Y[2*n+1]-Y[2*n-1]):
   s := sum(generalTerm, n=1..50);
   
   	
   s := 23.78942187
   
   
   The values are rather far apart. The second is not
   particularly close to -1, 0, or +1, values we might expect. Generate
   more data of your own by resetting _seed to your Social Security
   number (no hyphens please) and running
   for cntr from 1 to 10 do
   	X := {seq(r(), k=1..99)} union {0, 1}:
   	X := sort(convert(X, list)):
   	Y := {seq(r(), k=0..100)} union {0, 1}:
   	Y := sort(convert(Y, list)):
   	generalTerm := f(X[2*n],Y[2*n])
   		* (X[2*n+1]-X[2*n-1]) * (Y[2*n+1]-Y[2*n-1]):
   	s := sum(generalTerm, n=1..50):
   	print(s);
   od:
   
   	
   6.152782007
   	
   -6.223144712
   	
   -1.442629724
   	
   .5679381890
   	
   -.9000524038
   	
   1.302965578
   	
   -.7019236287
   	
   -4.458483906
   	
   -.1232440330
   	
   .9913757576
   
   
   
   Question 2 : Collect, either by running your procedure or gathering
   data from other groups, 100 values for Riemann sums of f over the
   unit square. Collect similar data for Riemann sums of g over the unit
   square and compare the two data sets. In particular, you may wish to
   pay attention to descriptive statistics such as mean, median, range,
   and standard deviation. (See ? stats,describe.)
   
   Question 3 : Are either of the double integrals of f or g over the unit
   square convergent? Argue for your conclusions.
   
   
   Extension
   Let C  be the unit cube in 3-space and define the function of three variables
   
   	
   
   
   Extend the variation of the Monte Carlo technique used above and calculate
   
   	
   
   
   Compare your results to that of a triple iterated integral.
   
   Report Requirements
   
   A minimal project report will include:
   
   	
    A written narrative supported by argumentation, data, and graphics for Questions
   		1 through 3.
   
   	
    In Question 2, also discuss how your data was collected and how you chose the
   		statistics to describe the data.