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

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 n^{3} in
R^{3}, after all!). Some sample problems are discussed here.

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=x^{2} 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])*

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.

We can also compute area via a double integral of 1dA as you learned in Chapter 15. To do a double integral, simply use oneint((2*x+3)-x^2,x=-1..3); 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.area:=int( int(1, y=x^2..(2*x+3)), x=-1..3); area:=32/3

Similarly, the total moment around the y-axis is the double integral of x over the region:x_moment:=int(int(y,y=x^2..2*x+3),x=-1..3); x_moment := 544/15

Then the center of gravity (here, due to the constant density, this is sometimes called they_moment:=int(int(x,y=x^2..2*x+3),x=-1..3); y_moment := 32/3

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 youcentroid:=[y_moment/area,x_moment/area]; centroid := [1, 17/5]

int(int(cos(x^2*y),y=x^2..2*x+3),x=-1..3);

This numerical approximation took 0.62 seconds on a fairly modern PC (running theevalf(int(int(cos(x^2*y),y=x^2..2*x+3),x=-1..3)); 3.326888539

Let's look at a region in space, inside the sphere
x^{2}+y^{2}+z^{2}=1 and "above" the paraboloid
z = x^{2} + y^{2}.

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

Sphere:=implicitplot3d(x^{2}+y^{2}+z^{2}=1,x=-1..1, y=-1..1, z=-1..1, color=blue, grid=[25,25,25], style=wireframe): Paraboloid:=plot3d(x^{2}+y^{2},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 5^{3}=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 R^{3} 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=x^{2}+y^{2} also satisfies
x^{2}+y^{2}+z^{2}=1. So we can substitute
z in for x^{2}+y^{2} in the second equation.
To find this value of z we enter:

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=xfsolve(z+z^2=1,z); -1.618033989, 0.6180339887

SphereTrimmed:=implicitplot3d(x^{2}+y^{2}+z^{2}=1,x=-1..1,y=-1..1,z=0.618..1, color=blue): ParaboloidTrimmed:=plot3d(x^{2}+y^{2},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);

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.

Note that the r in the integrand (the function that's integrated) comes from the Jacobian factor for cylindrical coordinates.bottom:=int(int(int(r,r=0..sqrt(z)),z=0..(0.6180339887)), theta=0..2*Pi); bottom := 0.5999908073

Above the intersection, the right-hand bound is given by the sphere r=sqrt(1-z

Hence we get the total volumetop:=int(int(int(r,r=0..sqrt(1-z^2)),z=0.6180339887..1), theta=0..2*Pi); top := 0.3999938717

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=rbottom+top; 0.9999846790

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.volume:=int(int(int(r,z=r^2..sqrt(1-r^2)),r=0..sqrt(0.6180339887)), theta=0..2*Pi); volume := 0.9999846791

Then a

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, andcyl_density:=subs({x=r*cos(theta),y=r*sin(theta)},density); cyl_density := z^{4}r^{2}cos(theta)^{2}+ 5 r^{2}sin(theta)^{2}

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.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

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.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

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