Information for Maple assignment 4

Each student will get individual data for the assignment.
Sketching regions in the plane
and computing double integrals 		Sketching regions in space
and computing triple integrals


Two dimensions

Integration in more than one dimension can pose difficulties which are new and unexpected compared to similar problems in one dimension. Certainly visualization of the regions can be difficult, and even numerical approximation of definite integrals can be computationally much lengthier (grid points go up as n3 in R3, after all!). Some sample problems are discussed here.



Drawing a picture of a region in the plane

Unfortunately and quite strangely, Maple does not have a built-in command to plot the region between two curves. There is a an option (filled=true) which will shade the region between a curve and the x-axis, but that is the extent of filling options. Hence we need to implement a slightly awkward work-around. We begin with something simple. The parabola y=x2 and the line y=2x+3 intersect at the points (-1,1) and (3,9). You can use solve(x2 = 2*x+3, x) (or similarly with fsolve) for these intersection points. The left-hand picture belowshows the boundary curves of the region, drawn with the command
plot({x^2,2*x+3},x=-2..4,thickness=2,scaling=constrained)

To plot the region between the curves, we will use plot3d to plot the plane z=0. We will not plot all of it, however. We will only plot the piece of z=0 bounded by the equations y=x2 and y=2x+3. We will set the orientation so that the initial view will an aerial view, making the region look like a two-dimensional plot.
plot3d(0, x=-1..3, y=x^2..(2*x+3), scaling=constrained, axes=normal, color=blue, style=patchnogrid, orientation=[-90,0])

If you play with orienation of the image, you will see what we are really plotting:

The choice between whether to use constrained or unconstrained scaling is up to the user. The choice should be made with the goal of effective visual communication in mind. You may wish to try both views in this and in other plotting situations, and choose which one seems best in the specific situation. This assignment commonly generates regions which are very tall and thin, which do not look good when scaling=constrained is used.

Computing the area of a region in the plane

We can compute the area of the region with either a single or a double integral. A single integral which gives the area is in the tradition of "Calc 1". 
 int((2*x+3)-x^2,x=-1..3);
                                 32/3
We can also compute area via a double integral of 1dA as you learned in Chapter 15. To do a double integral, simply use one int command inside another. 
area:=int( int(1, y=x^2..(2*x+3)), x=-1..3);
                                 area:=32/3
Notice how the ranges for the integrals are the same as the ranges when we plotted the region. This is not a coincidence.

Computing the center of gravity of a region in the plane

Suppose we think of the region enclosed by these two curves as a flat metal plate with constant density (what the text calls a lamina, see section 15.5). In fact, let's assume that the density is 1. Then the total moment of the lamina around the x-axis is given by taking the double integral of y over the region. 
 x_moment:=int(int(y,y=x^2..2*x+3),x=-1..3); 
                           x_moment := 544/15
Similarly, the total moment around the y-axis is the double integral of x over the region: 
 y_moment:=int(int(x,y=x^2..2*x+3),x=-1..3);
                          y_moment := 32/3
Then the center of gravity (here, due to the constant density, this is sometimes called the centroid) is given by "averaging" with respect to the total area: 
 centroid:=[y_moment/area,x_moment/area];
                        centroid := [1, 17/5]
We should then be able to balance the region on the tip of a finger at this point. You could try to "build" a region out of cardboard. Then you should be able to balance the region on the computed centroid or center of gravity. This may not be successful -- cutting things up and balancing them in the real world precisely may be more complicated than double integrals!

Computing a horrible double integral over a region in the plane

Out of curiousity, let's try to integrate some horrible function over this region, for example, cos(x2y). Let's try: 
int(int(cos(x^2*y),y=x^2..2*x+3),x=-1..3);
Maple gives up after finding the inside integral because it realizes that it can't find the antiderivative of these functions and so cannot use the Fundamental Theorem of Calculus. What if you really needed to know the value of this integral (or at least an approximate value)? You could do this: 
 evalf(int(int(cos(x^2*y),y=x^2..2*x+3),x=-1..3));
                                  3.326888539
This numerical approximation took 0.62 seconds on a fairly modern PC (running the Linux version of Maple). You might want to be cautious if you were studying fluid dynamics and needed approximate values of many integrals of even fairly simple functions over a six-dimensional region (space, momentum). Computational time might suddenly become important.

Three dimensions

Drawing a picture of a region in space

Let's look at a region in space, inside the sphere x2+y2+z2=1 and "above" the paraboloid z = x2 + y2.

I used these plotting commands after loading the plots package using the command with(plots):. 

Sphere:=implicitplot3d(x2+y2+z2=1,x=-1..1, y=-1..1, z=-1..1, color=blue, grid=[25,25,25], style=wireframe):
Paraboloid:=plot3d(x2+y2,x=-1..1,y=-1..1, color=green):
display({Sphere,Paraboloid}, axes=normal);

Take note of the options were specified for the plot of the sphere and you'll notice the grid option. Without going into too much detail on how the implicitplot3d command works, this grid option determines how smoothly to plot the surface. A higher grid setting means a smoother and more refined picture at the expense of longer computation time. The default is grid=[10,10,10] which results in a rather chunky sphere. For the curious, try grid=[5,5,5] to see how implicitplot3d sketches the surface. Using [50,50,50] is nice, of course, but computational time and storage go way up (roughly 53=125 times as much data needs to be computed and stored). The style=wireframe option allows us to see "inside" the sphere, so that the region of R3 under consideration can be inspected more closely.

Once you has a general idea of what the region looks like, you may want to take the extra step of "trimming" the plot. For example, we could plot only the portion of the sphere above the paraboloid, and only the portion of the paraboloid below the sphere. To do this, we'll need to know where the sphere and paraboloid intersect. The two surfaces intersect when the (x,y,z) triple which satisfies z=x2+y2 also satisfies x2+y2+z2=1. So we can substitute z in for x2+y2 in the second equation. To find this value of z we enter: 

 fsolve(z+z^2=1,z);
                          -1.618033989, 0.6180339887
Since the paraboloid opens "upwards" (positive values of z) the intersection will occur at (approximately!) 0.6180339887 and we don't need the other (negative) value of z. You can also confirm this by looking at the graph that was just displayed. Thus it will suffice to plot the sphere only in the range x=-1..1, y=-1..1, z=0.618..1. Plotting the paraboloid is a bit more complicated, as we must find the x and y ranges which keep the paraboloid beneath the sphere. The circle of intersection of the sphere and paraboloid is 0.6180339887=x2+y2, obtained by substituting the z-value into the equation of the paraboloid. Thus we may restrict the paraboloid to the region x=-sqrt(0.618)..sqrt(0.618), y=-sqrt(0.618-x^2)..sqrt(0.618-x^2). Putting all this information together we have: 
SphereTrimmed:=implicitplot3d(x2+y2+z2=1,x=-1..1,y=-1..1,z=0.618..1, color=blue):
ParaboloidTrimmed:=plot3d(x2+y2,x=-sqrt(0.618)..sqrt(0.618), y=-sqrt(0.618-x^2)..sqrt(0.618-x^2),color=green):
display3d({SphereTrimmed, ParaboloidTrimmed}, axes=normal, scaling=constrained);

Finding the volume of a region in space

The region is circularly symmetric and the z-axis is the axis of symmetry. The volume of the region can be computed using cylindrical coordinates. A cross-section is shown in the picture to the below. One can imagine our 3-d region as a solid of revolution of the radial cross-section below. The horizontal line is the intersection z=0.6180339887.

The roundness of the shapes suggests that cylindrical coordinates would be useful. The sphere is described by r2+z2=1 and the paraboloid described by z=r2. We may use the order of integration dr dz dθ or dz dr dθ. If we integrate with respect to r first, we will need to split the region into two parts: below the intersection z=0.618... and above the intesection. Below the intersection, the right-hand bound is given by the paraboloid r = sqrt(z). Thus the volume of the bottom is given by 
bottom:=int(int(int(r,r=0..sqrt(z)),z=0..(0.6180339887)), theta=0..2*Pi);
                            bottom := 0.5999908073
Note that the r in the integrand (the function that's integrated) comes from the Jacobian factor for cylindrical coordinates.

Above the intersection, the right-hand bound is given by the sphere r=sqrt(1-z2). Thus the volume of the bottom is given by 
top:=int(int(int(r,r=0..sqrt(1-z^2)),z=0.6180339887..1), theta=0..2*Pi);
                            top := 0.3999938717
Hence we get the total volume 
bottom+top;
                                 0.9999846790
If we used the order of integration dz dr dθ, then we only require one triple integral. The range for z is given below by the paraboloid z=r2 and above by the sphere z=sqrt(1-r^2). The range for r goes from r=0 (the z-axis) to whatever r solves z=r2 (equivalently r2+z2=1) when z=0.6180339887, that is, sqrt(0.6180339887). Thus the volume of the region is given by 
volume:=int(int(int(r,z=r^2..sqrt(1-r^2)),r=0..sqrt(0.6180339887)), theta=0..2*Pi);
                         volume := 0.9999846791
The difference in the tenth decimal place comes from compounding round-off errors in the two triple integrals when integrating with respect to r first.

Computing a mass given a density distribution

If the (mass) density at point (x,y,z) is given by (z4·x2)+5y2, the triple integral of the density over the whole region would be the total mass of the object. First, the density is stored to a variable.
density := z4 x2 + 5 y2
Then a substitution converts it into cylindrical coordinates. 
 cyl_density:=subs({x=r*cos(theta),y=r*sin(theta)},density);
                cyl_density := z4  r2  cos(theta)2  + 5 r2  sin(theta)2
We're still integrating in cylindrical coordinates, so the "extra" r (the Jacobian) is still needed. We can integrate with respect to r first and use two triple integrals, and Maple finds the mass as the sum of two triple iterated integrals.

bottom_mass:=int(int(int(cyl_density*r,r=0..sqrt(z)),z=0..(0.6180339887)), theta=0..2*Pi);
                            bottom_mass := 0.3128766261 

top_mass:=int(int(int(cyl_density*r,r=0..sqrt(1-z^2)),z=0.6180339887..1), theta=0..2*Pi);
                            top_mass := 0.2269499646 

bottom_mass+top_mass;
                                 0.5398265907
Thus the total mass is about 0.5398 units. If we integrate with respect to z first, we need only one triple integral to get the mass. 
mass := int(int(int(cyl_density*r,z=r^2..sqrt(1-r^2)),r=0..sqrt(0.6180339887)), theta=0..2*Pi);
                         mass := 0.5398265909
This confirms that the total mass is about 0.5398 units. The discrepancy in the tenth decimal place is again most likely due to compounding round-off errors in the first estimation.

Maintained by greenfie@math.rutgers.edu and last modified 6/14/2008 by Andrew Baxter.