Information for Maple assignment 4
Each student will get individual data by
e-mail for the assignment or the data will be posted on a web page.
Two dimensions
Integration in dimensions >1 can pose difficulties which are new
and unexpected compared to similar problems in 1 dimension. Certainly
visualization of the regions can be difficult, and even
numerical approximation of definite integrals can be computationally
much more difficult (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
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 solve or fsolve for
this information.) The first picture to the right shows the boundary
curves of the region, drawn with the command
plot({x^2,2*x+3},x=-2..4,thickness=2,scaling=constrained)
The package plots has a command that
allows the region itself to be viewed.
plot3d(0,y=x^2..(2*x+3),x=-1..3,scaling=constrained, axes=normal,color=blue,style=patchnogrid,orientation=[0,179])
The various options show the picture displayed second to the right.
Please note that the x boundaries in the two commands are
different. Also please note that the option scaling=constrained is used in both
cases. The choice between the constrained and unconstrained view is up
to the user. The choice should be made with the
goal of effective visual communication. You may wish to try
both views in this and in other plotting situations, and choose which
one seems best in the specific situation.
|
|
|
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. The single integral which gives the area is in the tradition
of "Calc 1". The double integral is the sum of 1 dA's over the
region, and is more the style of several variable calculus.> int((2*x+3)-x^2,x=-1..3);
32/3
> area:=int(int(1,y=x^2..(2*x+3)),x=-1..3);
area := 32/3
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
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);
544
x_moment := ---
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
You can even 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);
bytes used=4002008, alloc=3407248, time=0.13
bytes used=8002304, alloc=5307444, time=0.28
bytes used=12002576, alloc=5897160, time=0.44
3
/ 4 2
| sin(x ) - sin(x (2 x + 3))
| - --------------------------- dx
| 2
/ x
-1
Maple gives up (quite rapidly,
actually!) because it realizes that it can't find the antiderivative
of these functions and so it can't use the Fundamental Theorem of
Calculus. The various "messages" shown signal the user about the storage
and computation time. What if you really needed to know a value
(at least, an approximate value) of this integral? 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 appropriate package
using the command with(plots):.
>A:=plot3d(x2+y2,x=-1..1,y=-1..1,axes=normal,color=green):
>B:=implicitplot3d(x2+y2+z2=1,x=-1..1,y=-1..1,z=-1..1,
color=blue,axes=normal,grid=[25,25,25],style=wireframe):
>display3d({A,B});
Several options were specified for the plot of the sphere. grid=[25,25,25] requests a denser sampling
than usual. The default is [10,10,10] which results in a spherical
object that isn't entirely "convincing" to me. For amusement, you
might try [5,5,5] which shows 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.
|
 |
Finding the volume of a region in space
The two surfaces intersect when the (x,y,z) triple which satisfies
z=x2+y2 also satisfies
x2+y2+z2=1. So we can "plug"
z in for x2+y2 in the second equation.
To find this value of z we write:> fsolve(z^2+z=1,z);
-1.618033989, 0.6180339887
Since the paraboloid opens "up" (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.
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
left. There is one integral as z goes from 0 to 0.6180339887 with the
outer boundary equal to the paraboloid: z=x2+y2,
or, in cylindrical coordinates, z=r2 which is r=sqrt(z). The
second integral has z going from 0.6180339887 to 1 with the outer
boundary determined by the sphere: x2+y2+z2=1, or, in cylindrical
coordinates, r2+z2=1 which is r=sqrt(1-z2)). The
r inside (the integrand or function
that's integrated) comes from the Jacobian factor
for cylindrical coordinates. >part_1:=int(int(int(r,r=0..sqrt(z)),theta=0..2*Pi),z=0..(0.6180339887));
part_1 := 0.5999908073
> part_2:=int(int(int(r,r=0..sqrt(1-z^2)),theta=0..2*Pi),z=(0.6180339887)..1);
part_2 := 0.3999938717
> part_1+part_2;
0.9999846790
This total is the volume of the region inside the sphere and above the paraboloid.
|
|
Computing a mass given a density distribution
If the region had a mass density given by
(z4·x2)+5y2, the triple integral of
the mass density over the whole region would be the total mass of the
object. First, the density is defined, and then a substitution converts it into cylindrical
coordinates. The "extra" r is still needed (the Jacobian) and
Maple finds the mass as the sum of two
triple iterated integrals.> rho:=z^4*x^2+5*y^2;
4 2 2
rho := z x + 5 y
> cyl_rho:=subs({x=r*cos(theta),y=r*sin(theta)},rho);
4 2 2 2 2
cyl_rho := z r cos(theta) + 5 r sin(theta)
> part_1:=int(int(int(cyl_rho*r,r=0..sqrt(z)),theta=0..2*Pi),z=0..(0.6180339887));
part_1 := 0.3128766268
> part_2:=int(int(int(cyl_rho*r,r=0..sqrt(1-z^2)),theta=0..2*Pi),z=(0.6180339887)..1);
part_2 := 0.2269499638
> part_1+part_2;
0.5398265906
Maintained by
greenfie@math.rutgers.edu and last modified 11/8/2006.