Multiple Integrals 

In this worksheet we will explore computing multiple integrals in rectangular, cylindrical, and spherical coordinates. In addition, there are a few examples using Jacobians toward the end of the file. For the most part we will stick to using the "solve" command to help us determine bounds and then used nested copies of "int" to compute iterated integrals. 

Maple worksheet: multiple_integrals.mw 

>	restart; with(plots): with(VectorCalculus):

 

Example: We can parametrize any surface of the form: using two parameters (say and ) as follows: , , and . We can then plot such a surface using either in the "old" way or the "new" (parametric) way.  

Consider the function plotted over the square . The top plot3d command plots the graph of the "old way" and the bottom plot3d command plots the graph thought as a parametric surface [note since we used the rectangular coordinate system to parameterize our surface, the plots look identical]. 

>	f := (x,y) -> 4 - x^2/8 - y^2/8: 'f(x,y)' = f(x,y); plot3d(f(x,y),x=-5..5,y=-5..5,axes=boxed); plot3d([s,t,f(s,t)],s=-5..5,t=-5..5,axes=boxed);

 

We can compute the volume under this paraboloid using the iterated integral   

>	int(int(f(x,y),y=-5..5),x=-5..5);

 

(1)

 

Maple can handle plots defined over "x-simple" regions (regions whose left & right sides are graphs of functions) or "y-simple" regions (regions whose top & bottoms are graphs of functions). For example, let's graph over the region (inside a circle centered at the origin with radius 4). This region is both x and y-simple. Let's consider it as x-simple first. We find bounds for in terms of and then constant bounds.  

Although I will use Maple to solve for these bounds, this is certainly overkill. We could easy find them "by hand" or by merely looking at it and writing down the "obvious" answers. 

>	xBounds := solve(x^2+y^2=16,x);

 

(2)

 

>	solve({xBounds},y);

 

(3)

 

So we get the bounds: and . 

>	plot3d(f(x,y),x=-sqrt(16-y^2)..sqrt(16-y^2),y=-4..4,scaling=constrained,axes=boxed,orientation=[45,55,0]);

 

Now let's compute the volume under using the double integral and these bounds. 

>	Int(Int(f(x,y),x=-sqrt(16-y^2)..sqrt(16-y^2)),y=-4..4) = int(int(f(x,y),x=-sqrt(16-y^2)..sqrt(16-y^2)),y=-4..4);

 

(4)

 

Of course it would make far more sense to convert to polar coordinates and compute the integral that way... 

>	Int(Int((4-r^2/8)*r,r=0..4),theta=0..2*Pi) = int(int(f(r*cos(theta),r*sin(theta))*r,r=0..4),theta=0..2*Pi);

 

(5)

 

Example: Find the average value of over the region inside the sphere centered at the origin of radius 2 lying in the first octant. Note that the average value of over a region E is -- that is -- the integral of over E divided by the volume of E. 

>	f := (x,y,z) -> x*y*z: 'f(x,y,z)' = f(x,y,z);

 

(6)

 

To plot the region of integral let's plot the surfaces which make up the sides of region. For the part of the ball of radius 2 lying in the first octant, I need the sphere of radius 2 (this makes up the top, front, and right side), part of the -plane (the bottom), part of the -plane (the back), and part of the -plane (the left side). 

>

sPlot := plot3d(sqrt(4-x^2-y^2),x=0..2,y=0..sqrt(4-x^2),color=blue): xyPlot := plot3d(0,x=0..2,y=0..sqrt(4-x^2),color=green): yzPlot := plot3d([0,y,z],y=0..sqrt(4-z^2),z=0..2,color=red): xzPlot := plot3d([x,0,z],x=0..2,z=0..sqrt(4-x^2),color=gray): display({sPlot,xyPlot,yzPlot,xzPlot},scaling=constrained,axes=boxed,orientation=[-115,-40,35]);

 

>	vol := int(int(int(1       , z=0..sqrt(4-x^2-y^2)), y=0..sqrt(4-x^2)), x=0..2);

 

(7)

 

>	fint := int(int(int(f(x,y,z), z=0..sqrt(4-x^2-y^2)), y=0..sqrt(4-x^2)), x=0..2);

 

(8)

 

>	"Average Value" = fint / vol;

 

(9)

 

>	Int(Int(Int(f(x,y,z), z=0..sqrt(4-x^2-y^2)), y=0..sqrt(4-x^2)), x=0..2) / Int(Int(Int(1       , z=0..sqrt(4-x^2-y^2)), y=0..sqrt(4-x^2)), x=0..2) = fint / vol;

 

(10)

 

Note: To explain the last "weird" command. "Int" is called an inert command -- Maple leaves the integral inert -- it does not actually compute the integral. It just prints it out. In general if I want to display what I'm computing, I like to type something like "Inert commands = normal commands" so that Maple prints what I'm computing on the left hand side of the equation and actually computes the answers for the right hand side. For example: 

>	Int(sin(x),x=0..Pi) = int(sin(x),x=0..Pi);

 

(11)

 

Example: Consider the sphere centered at the origin with radius 5. In rectangular coordinates it is described by , but in spherical coordinates we have the much simpler equation "" (rho is equal to 5). Using these coordinates, we see that:, , and  where      and   . Let's plot our sphere using a 3D parametric plot. 

>	plot3d([5*cos(theta)*sin(phi),5*sin(theta)*sin(phi),5*cos(phi)],theta=0..2*Pi,phi=0..Pi,axes=boxed);

 

This looks a lot nicer than what we would get with "implicitplot3d" or plotting the upper and lower halves together using plot3d... 

>	implicitplot3d(x^2+y^2+z^2=5^2,x=-5..5,y=-5..5,z=-5..5,axes=boxed); display({ plot3d( sqrt(5^2-x^2-y^2),x=-5..5,y=-sqrt(5^2-x^2)..sqrt(5^2-x^2)),          plot3d(-sqrt(5^2-x^2-y^2),x=-5..5,y=-sqrt(5^2-x^2)..sqrt(5^2-x^2)) }, axes=boxed, orientation=[-55,65,15]);

 

Now let's compute the volume of the sphere first using rectangular coordinates (whose bounds come from the pair of plot3d plots) and then spherical coordinates (whose bounds come from the parametric plot). 

>	int(int(int(1,z=-sqrt(5^2-x^2-y^2)..sqrt(5^2-x^2-y^2)),y=-sqrt(5^2-x^2)..sqrt(5^2-x^2)),x=-5..5);

 

(12)

 

>	int(int(int(1*rho^2*sin(phi),rho=0..5),phi=0..Pi),theta=0..2*Pi);

 

(13)

 

Example: Let  be the region in bounded by and . 

First, let's make a plot of  (Well, actually, we'll plot the outer surface of the solid region ). 

Notice that these two surfaces intersect when . So we should plot over the disk of radius 2 centered at the origin in the -plane. The bottom of this region is given by and the top by . These semicircular arcs range from to . So the bounds for x and y in our plot3d commands should be "y=-sqrt(4-x^2)..sqrt(4-x^2),x=-2..2". The top of the region is given by the horizontal plane "" and the bottom is given by the cone "". 

>	top := plot3d(2,y=-sqrt(4-x^2)..sqrt(4-x^2),x=-2..2):

 

>	bottom := plot3d(sqrt(x^2+y^2),y=-sqrt(4-x^2)..sqrt(4-x^2),x=-2..2):

 

>	display({top,bottom});

 

Alternatively, we could plot these surfaces parametrically. The most natural coordinates to use are "cylindrical coordinates". In this case is still , but   becomes .  We have x and y ranging over the interior of the disk centered at the origin of radius 2. In cylindrical coordinates, translates to and . Thus the top surface translates to , , . The bottom surface (the cone) translates to , , . 

>	top := plot3d([r*cos(theta),r*sin(theta),2],r=0..2,theta=0..2*Pi): bottom := plot3d([r*cos(theta),r*sin(theta),r],r=0..2,theta=0..2*Pi): display({top,bottom});

 

Notice that translating to cylindrical coordinates yields a nicer looking plot (because we're using the "proper" coordinate system). 

Next, let's compute the triple integral . We need to write it as an iterated integral. 

First, we will compute this integral using rectangular coordinates (and do it in all 6 orders of integration). Of course, you we should get the SAME answer each time so this is kind of a silly thing to do. [WARNING: Maple kind of "chokes" on some of these integrals. We have to add an "evalf" command to get a "real" answer.] 

For the first two integrals, we focus on the top and bottom surfaces. The lower surface is   and upper surface is . So these provide our bounds for the first two integrals. The and bounds are the same as those in the plot3d commands. 

The next two integrals require bounds first. So we need to find the back and front of the region. Solving   for we get . To get the and bounds we need to consider what the region "squished" onto the -plane would look like. Looking the plots, we can see this is just the slice of the region when cut through by the -plane. Setting we get which is . Intersecting this with we get -limits of . So if comes first, our bounds are and . If comes first we need to integrate from the left branch of the absolute value function to the right branch. This gives us the bounds and . 

The last two integrals are just like the middle two integrals (we just need to interchange x's and y's). 

>	plot({abs(y),2},y=-2..2,labels=[y,z],color=[red,red],title="The projection of E onto the yz-plane");

 

>	f := (x,y,z) -> x^2+y*z: 'f(x,y,z)' = f(x,y,z);

 

(14)

 

>	Int(Int(Int(f(x,y,z),z=sqrt(x^2+y^2)..2),y=-sqrt(4-x^2)..sqrt(4-x^2)),x=-2..2) = evalf(int(int(int(f(x,y,z),z=sqrt(x^2+y^2)..2),y=-sqrt(4-x^2)..sqrt(4-x^2)),x=-2..2));

 

(15)

 

>	Int(Int(Int(f(x,y,z),z=sqrt(x^2+y^2)..2),x=-sqrt(4-y^2)..sqrt(4-y^2)),y=-2..2) = evalf(int(int(int(f(x,y,z),z=sqrt(x^2+y^2)..2),x=-sqrt(4-y^2)..sqrt(4-y^2)),y=-2..2));

 

(16)

 

>	Int(Int(Int(f(x,y,z),x=-sqrt(z^2-y^2)..sqrt(z^2-y^2)),y=-z..z),z=0..2) = int(int(int(f(x,y,z),x=-sqrt(z^2-y^2)..sqrt(z^2-y^2)),y=-z..z),z=0..2);

 

(17)

 

>	Int(Int(Int(f(x,y,z),x=-sqrt(z^2-y^2)..sqrt(z^2-y^2)),z=abs(y)..2),y=-2..2) = evalf(int(int(int(f(x,y,z),x=-sqrt(z^2-y^2)..sqrt(z^2-y^2)),z=abs(y)..2),y=-2..2));

 

(18)

 

>	Int(Int(Int(f(x,y,z),y=-sqrt(z^2-x^2)..sqrt(z^2-x^2)),x=-z..z),z=0..2) = int(int(int(f(x,y,z),y=-sqrt(z^2-x^2)..sqrt(z^2-x^2)),x=-z..z),z=0..2);

 

(19)

 

>	Int(Int(Int(f(x,y,z),y=-sqrt(z^2-x^2)..sqrt(z^2-x^2)),z=abs(x)..2),x=-2..2) = evalf(int(int(int(f(x,y,z),y=-sqrt(z^2-x^2)..sqrt(z^2-x^2)),z=abs(x)..2),x=-2..2));

 

(20)

 

So  

Of course this is kind of a silly way to compute this triple integral. Cylindrical coordinates should make this much easier to compute. Let's convert the integral to cylindrical coordinates and compute the answer. First, let's translate our function... 

>	"f in cylindrical coordinates" = f(r*cos(theta),r*sin(theta),z);

 

(21)

 

Next, the lower bound for is and the upper bound is . Then in the -plane our bounds are the same as when we plotted the region using parameterized surfaces: and . [Don't forget the Jacobian "r"!!] 

>	Int(Int(Int(f(r*cos(theta),r*sin(theta),z)*r,z=r..2),r=0..2),theta=0..2*Pi)= int(int(int(f(r*cos(theta),r*sin(theta),z)*r,z=r..2),r=0..2),theta=0..2*Pi);

 

(22)

 

Next, let's convert the integral to spherical coordinates and compute the answer (yet again). It's easy to see that we're are integrating "all the way around" so that "" (just like cylindrical coordinates). Also, ranges from 0 to the side of the cone. Converting the cone  

to spherical coordinates, we get . So we must have which tells us that 's upper bound is . Finally, notice that should range from 0 until we hit the plane . Converting this to spherical coordinates, we get. So the upper bound for is . [Don't forget the Jacobian ""!!] 

>	"f in spherical coordinates" = f(rho*cos(theta)*sin(phi),rho*sin(theta)*sin(phi),rho*cos(phi));

 

(23)

 

>	Int(Int(Int(f(rho*cos(theta)*sin(phi),rho*sin(theta)*sin(phi),rho*cos(phi))*rho^2*sin(phi),rho=0..2*sec(phi)),theta=0..2*Pi),phi=0..Pi/4)= int(int(int(f(rho*cos(theta)*sin(phi),rho*sin(theta)*sin(phi),rho*cos(phi))*rho^2*sin(phi),rho=0..2*sec(phi)),theta=0..2*Pi),phi=0..Pi/4);

 

(24)

 

Example: Let's plot the region of integration of the integral: and then use that to rewrite the integral in the other 5 orders of intergration. The bounds on the integrals tell us that we are integrating over the solid region bounded by  ,  ,   ,  and  . [It turns out that the other bounds (the upper bounds on the outer two integrals): and are extraneous. If we plot them along with the bounds listed above, they will just make it harder to "see" what's going on.]  Let's simultaneously plot all of the relevant surfaces and then make the proper restrictions to get a plot of our "mystery" solid.  

First, we'll create plots for the surfaces and using a standard plot3d command (since these are graphs of functions of x and y). We'll restrict our inputs as follows: (this seems like a reasonable first guess given the bounds for x and y in the integral above). Notice that graphs of functions of x and y should give us a "top" and "bottom". 

>	bottom := plot3d(0,x=0..1,y=0..1,color=blue): top := plot3d(1-y,x=0..1,y=0..1,color=green):

 

Next, if we wish to find a "front" and a "back" we need surfaces be of the form: (if we wanted to plot "left and right sides", we should have surfaces of the form: ). Let's plot (the back) and (the front -- we got this by solving for ). Again we'll restrict our inputs y and z to the square: because it seems like a reasonable place to start. 

Now this time we can't just use a standard plot3d command (it handles type graphs). Instead these new plots must be done "parametrically". We let , , and . So first, and then . 

>	back := plot3d([0,y,z],y=0..1,z=0..1,color=red): front := plot3d([y^2,y,z],y=0..1,z=0..1,color=yellow):

 

>	display({bottom,top,back,front},axes=normal,labels=[x,y,z],viewpoint=circleleft);

 

From this we can see how to "cut off" irrelevant parts of our surfaces. 

>	bottom := plot3d(0,x=0..1,y=sqrt(x)..1,color=blue): top := plot3d(1-y,x=0..1,y=sqrt(x)..1,color=green): back := plot3d([0,y,z],y=0..1,z=0..1-y,color=red): front := plot3d([y^2,y,z],y=0..1,z=0..1-y,color=yellow): display({bottom,top,back,front},axes=normal,labels=[x,y,z],viewpoint=circleleft);

 

Using this graph, we can better "see" how to write down the remaining 5 iterated integrals. For the rest of the example, let's assume the function is and then compute all 6 iterated integrals. [If our answers are right, we should get the same result each time.] 

First, "squash out" the -direction and see what region in the -plane we are dealing with (go back to out "cut-down" plot and look at it from "above"). 

>	display({bottom,top,back,front},axes=normal,labels=[x,y,z],orientation=[270,0,0]);

 

So the y bounds are and then .  Or in the other order, the x bounds are and . 

>	Int(Int(Int(x^2+y*z,z=0..1-y),y=sqrt(x)..1),x=0..1) = int(int(int(x^2+y*z,z=0..1-y),y=sqrt(x)..1),x=0..1); Int(Int(Int(x^2+y*z,z=0..1-y),x=0..y^2),y=0..1) = int(int(int(x^2+y*z,z=0..1-y),x=0..y^2),y=0..1);

 

 

	(25)

 

In the next two integrals we will place first, so let's see what region of the-plane we are dealing with when the 's are "squashed out". 

>	display({bottom,top,back,front},axes=normal,labels=[x,y,z],orientation=[-90,0,90]);

 

So after integrating from to (eliminating ) we need to let range from 0 to and then to range from 0 to 1. Interchanging y and z yields similar bounds. 

>	Int(Int(Int(x^2+y*z,x=0..y^2),z=0..1-y),y=0..1) = int(int(int(x^2+y*z,x=0..y^2),z=0..1-y),y=0..1); Int(Int(Int(x^2+y*z,x=0..y^2),z=0..1-y),y=0..1) = int(int(int(x^2+y*z,x=0..y^2),y=0..1-z),z=0..1);

 

 

	(26)

 

In the final two integrals we will place first, so let's see what region of the-plane we are dealing with when the 's are "squashed out". 

>	display({bottom,top,back,front},axes=normal,labels=[x,y,z],orientation=[180,90,-90]);

 

ranges from the yellow surface to green surface ( to ). In the -plane, x and z range from 0 to the curve determined by where the yellow and green surfaces intersect and . So and thus and . 

>	Int(Int(Int(x^2+y*z,y=sqrt(x)..1-z),z=0..1-sqrt(x)),x=0..1) = int(int(int(x^2+y*z,y=sqrt(x)..1-z),z=0..1-sqrt(x)),x=0..1); Int(Int(Int(x^2+y*z,y=sqrt(x)..1-z),x=0..(1-z)^2),z=0..1) = int(int(int(x^2+y*z,y=sqrt(x)..1-z),x=0..(1-z)^2),z=0..1);

 

 

	(27)

 

Example: Consider the solid region E bounded by , , , , and and let  .    

>	f := (x,y,z) -> x+y^2+z^3: 'f(x,y,z)' = f(x,y,z);

 

(28)

 

First, we will plot the 5 surfaces (in the same plot window) which define the bounded region E. The top and bottom and can be plotted using the standard plot3d command. The other surfaces have to be treated as parametric surfaces. For example, is handled using the parametrization , , and . In so in the plot3d command we have "[x,1-x,z]". 

>

top := plot3d(1-x^2,x=0..1,y=0..1,color=blue): bottom := plot3d(0,x=0..1,y=0..1,color=green): front := plot3d([x,1-x,z],x=0..1,z=0..1,color=red): back := plot3d([0,y,z],y=0..1,z=0..1,color=yellow): side := plot3d([x,0,z],x=0..1,z=0..1,color=grey): display({top,bottom,front,back,side}, viewpoint=circleleft);

 

Now we get rid of the unnecessary parts of the boundary surfaces. We can find these bounds by intersecting pairs of surfaces. For example, the blue surface is cut by the red surface. So provides us with a bound for the blue surface.   

>

top := plot3d(1-x^2,x=0..1,y=0..1-x,color=blue): bottom := plot3d(0,x=0..1,y=0..1-x,color=green): front := plot3d([x,1-x,z],x=0..1,z=0..1-x^2,color=red): back := plot3d([0,y,z],y=0..1,z=0..1,color=yellow): side := plot3d([x,0,z],x=0..1,z=0..1-x^2,color=grey): display({top,bottom,front,back,side},viewpoint=circleleft);

 

Let's compute the integral  in all 6 orders of integration. We'll start with the integrals with dz coming first. 

The z bounds are and . Next, we'll look at the region of integration after "squishing out" the z direction. 

>	display({bottom,top,back,front,side},axes=normal,labels=[x,y,z],orientation=[270,0,0]);

 

So x and y are bounded by 0 below and above. 

>	Int(Int(Int(f(x,y,z),z=0..1-x^2),x=0..1-y),y=0..1) = int(int(int(f(x,y,z),z=0..1-x^2),x=0..1-y),y=0..1); Int(Int(Int(f(x,y,z),z=0..1-x^2),y=0..1-x),x=0..1) = int(int(int(f(x,y,z),z=0..1-x^2),y=0..1-x),x=0..1);

 

 

	(29)

 

In the final two integrals we will place first, so let's see what region of the-plane we are dealing with when the 's are "squashed out". 

>	display({bottom,top,back,front,side},axes=normal,labels=[x,y,z],orientation=[180,90,-90]);

 

The blue surface provides an upper bound for x and z that is . 

>	Int(Int(Int(f(x,y,z),y=0..1-x),x=0..sqrt(1-z)),z=0..1) = int(int(int(f(x,y,z),y=0..1-x),x=0..sqrt(1-z)),z=0..1); Int(Int(Int(f(x,y,z),y=0..1-x),z=0..1-x^2),x=0..1) = int(int(int(f(x,y,z),y=0..1-x),z=0..1-x^2),x=0..1);

 

 

	(30)

 

In the next two integrals we will place first, so let's see what region of the-plane we are dealing with when the 's are "squashed out". 

>	display({bottom,top,back,front,side},axes=normal,labels=[x,y,z],orientation=[270,180,270]);

 

Notice that even though y and z are ranging over the square , it is broken into 2 pieces. So we will need to break our integral into 2 parts. 

The yellow surface provides a lower bound for both pieces. The upper bound is given by (blue) and  (red). To find bounds for y and z we need to intersect these 2 surfaces: and so and thus and also . 

>

Int(Int(Int(f(x,y,z),x=0..sqrt(1-z)),y=0..1-sqrt(1-z)),z=0..1)+Int(Int(Int(f(x,y,z),x=0..1-y),y=1-sqrt(1-z)..1),z=0..1) = int(int(int(f(x,y,z),x=0..sqrt(1-z)),y=0..1-sqrt(1-z)),z=0..1)+int(int(int(f(x,y,z),x=0..1-y),y=1-sqrt(1-z)..1),z=0..1); Int(Int(Int(f(x,y,z),x=0..sqrt(1-z)),z=1-(1-y)^2..1),y=0..1)+Int(Int(Int(f(x,y,z),x=0..1-y),z=0..1-(1-y)^2),y=0..1) = int(int(int(f(x,y,z),x=0..sqrt(1-z)),z=1-(1-y)^2..1),y=0..1)+int(int(int(f(x,y,z),x=0..1-y),z=0..1-(1-y)^2),y=0..1);

 

 

	(31)

 

Example:  Let and consider the region which is bounded by and . 

>	f := (x,y,z) -> x^2+y*z: 'f(x,y,z)' = f(x,y,z);

 

(32)

 

To plot the region of integral I will plot the surfaces which make up the sides of region. This consists of the xy, xz, and yz-planes and two planes which form a "roof" for .  

Note: Be careful: regular "plot3d(function, ranges);" doesn't really pay attention to the parameter names you use like "x=...." and "y=...."  It will use the first range as a range of "x's" and the second range as a range of "y's". So, for example, "plot3d(1, y=0..2, x=0..1);  is not the same as plot3d(1, x=0..1, y=0..2); The second command has x and y in the "correct" order. 

A few notes on where the bounds for various variables come from... 

The gray side's boundary is determined by line of intersection of the green and blue planes projected onto the xz-plane. The blue plane's equation is y=2-2z and the green's is z=1-x. So in the xz-plane they cross at z=1-x. Thus the bounds for the gray side should change to "z=0..1-x" 

The green and blue plane cross where y=2-2z and z=1-x both hold. To find green bounds we need to project this intersection into the xy-plane (the green plane's bounds are given in terms of x and y). Thus subbing z=1-x into y=2-2z, we get y=2-2(1-x) so that y=2x. So the bounds for y (for the green plane) should start at y=0 and go up to y=2x. 

The blue plane is drawn using the parameters x and z. So we should eliminate y from the equations z=1-x and y=2-2z. But this is already done in the first equation: z=1-x. So z should range from 0 to 1-x. [This is the same as the gray plane.] 

>

xyPlot := plot3d(0,x=0..1,y=0..2,color=pink): yzPlot := plot3d([0,y,z],y=0..2-2*z,z=0..1,color=red): xzPlot := plot3d([x,0,z],x=0..1,z=0..1-x,color=gray): top1Plot := plot3d([x,2-2*z,z],x=0..1,z=0..1-x,color=blue): top2Plot := plot3d(1-x,x=0..1,y=0..2*x,color=green): display({xyPlot,yzPlot,xzPlot,top1Plot,top2Plot},scaling=constrained,axes=boxed,viewpoint=circleleft);

 

Next, let's compute this integral using all 6 orders of integration. We can use the bounds found while fixing the graph above to help with the integrals below. Notice that the orders "dz dy dx" and "dz dx dy" require us to split the integral into 2 pieces (1 with a "blue roof" and 1 with a "green roof"). The first two integrals go from bottom (z=0) to top. However, the top has two pieces. The green piece: z=1-x and the blue piece: z=1-y/2. After z is integrated out, the piece of the xy-plane the solid sits over must be divided into two pieces (two triangles). These 

triangles are separated by the line y=2x (the intersection of the planes projected down). 

>	int(int(int(f(x,y,z),z=0..1-x), y=0..2*x), x=0..1) + int(int(int(f(x,y,z),z=0..1-y/2), y=2*x..2), x=0..1);

 

(33)

 

>	int(int(int(f(x,y,z),z=0..1-x), x=y/2..1), y=0..2) + int(int(int(f(x,y,z),z=0..1-y/2), x=0..y/2), y=0..2);

 

(34)

 

>	int(int(int(f(x,y,z),y=0..2-2*z), z=0..1-x), x=0..1);

 

(35)

 

>	int(int(int(f(x,y,z),y=0..2-2*z), x=0..1-z), z=0..1);

 

(36)

 

>	int(int(int(f(x,y,z),x=0..1-z), z=0..1-y/2), y=0..2);

 

(37)

 

>	int(int(int(f(x,y,z),x=0..1-z), y=0..2-2*z), z=0..1);

 

(38)

 

Example: Let's plot the surfaces which bound the region of integration where is bounded by and and .  

By intersecting and we get . So the paraboloid should be plotted over a region inside the circle of radius 2. Noting that , we have the bounds and then . 

>	top := plot3d(4-x^2-y^2,x=-2..2,y=0..sqrt(4-x^2),color=red): bottom := plot3d(0,x=-2..2,y=0..sqrt(4-x^2),color=green): side := plot3d([x,0,z],x=-2..2,z=0..4-x^2,color=blue): display({bottom,top,side},scaling=constrained,axes=boxed,viewpoint=circleleft);

 

Let's compute the integral by using the 6 different orders of integration. 

>	Int(Int(Int(x^2+y^2, z=0..4-x^2-y^2), y=0..sqrt(4-x^2)), x=-2..2) = int(int(int(x^2+y^2, z=0..4-x^2-y^2), y=0..sqrt(4-x^2)), x=-2..2);

 

>	Int(Int(Int(x^2+y^2, z=0..4-x^2-y^2), x=-sqrt(4-y^2)..sqrt(4-y^2)), y=0..2) = int(int(int(x^2+y^2, z=0..4-x^2-y^2), x=-sqrt(4-y^2)..sqrt(4-y^2)), y=0..2);

 

>	Int(Int(Int(x^2+y^2, y=0..sqrt(4-x^2-z)), z=0..4-x^2), x=-2..2) = int(int(int(x^2+y^2, y=0..sqrt(4-x^2-z)), z=0..4-x^2), x=-2..2);

 

>	Int(Int(Int(x^2+y^2, y=0..sqrt(4-x^2-z)), x=-sqrt(4-z)..sqrt(4-z)), z=0..4) = int(int(int(x^2+y^2, y=0..sqrt(4-x^2-z)), x=-sqrt(4-z)..sqrt(4-z)), z=0..4);

 

>	Int(Int(Int(x^2+y^2, x=-sqrt(4-y^2-z)..sqrt(4-y^2-z)), z=0..4-y^2), y=0..2) = int(int(int(x^2+y^2, x=-sqrt(4-y^2-z)..sqrt(4-y^2-z)), z=0..4-y^2), y=0..2);

 

>	Int(Int(Int(x^2+y^2, x=-sqrt(4-y^2-z)..sqrt(4-y^2-z)), y=0..sqrt(4-z)), z=0..4) = int(int(int(x^2+y^2, x=-sqrt(4-y^2-z)..sqrt(4-y^2-z)), y=0..sqrt(4-z)), z=0..4);

 

 

 

 

 

 

	(39)

 

This problem is more naturally set in terms of cylindrical coordinates. In cylindrical coordinates, we have , , and . Also, our function simplifies to (and don't forget the Jacobian ""). Let's plot the region again using parametric descriptions coming from cylindrical coordinates and then recompute the integral. Notice that by choosing the "correct" coordinate system, our plot looks a bit nicer and Maple doesn't have to work as hard computing the integral. 

>

top := plot3d([r*cos(theta),r*sin(theta),4-r^2],r=0..2,theta=0..Pi,color=red): bottom := plot3d([r*cos(theta),r*sin(theta),0],r=0..2,theta=0..Pi,color=green): side := plot3d([r*cos(theta),0,4-r^2],r=0..2,theta=0..Pi,color=blue): display({bottom,top,side},scaling=constrained,axes=boxed,viewpoint=circleleft); Int(Int(Int(r^2*r,z=0..4-r^2),r=0..2),theta=0..Pi) = int(int(int(r^2*r,z=0..4-r^2),r=0..2),theta=0..Pi);

 

 

	(40)

 

Example: Let and consider the region which is bounded by and . 

>	f := (x,y,z) -> x^2+y*exp(z): 'f(x,y,z)' = f(x,y,z);

 

(41)

 

To plot the region of integral we will plot the surfaces which make up the sides of region. This consists of the xy and xz-planes as well asand two more surfaces which form the "roof" of . As before, intersecting surfaces helps us find bounds. 

>	solve({z=x^2,y=0}); solve({z=x^2,y=3-3*z});

 

 

	(42)

 

>

xyPlot := plot3d(0,x=0..1,y=0..3,color=pink): yzPlot := plot3d([1,y,z],y=0..3-3*z,z=0..1,color=red): xzPlot := plot3d([x,0,z],x=0..1,z=0..x^2,color=gray): top1Plot := plot3d([x,3-3*z,z],x=0..1,z=0..x^2,color=blue): top2Plot := plot3d(x^2,x=0..1,y=0..3-3*x^2,color=green): display({xyPlot,yzPlot,xzPlot,top1Plot,top2Plot},scaling=constrained,axes=boxed,viewpoint=circleleft);

 

Now let's compute this integral using all 6 orders of integration.  

Note: For the orders of integration "dz dy dx" and "dz dx dy", the we will need to split the integral into 2 pieces (1 with a "blue roof" and 1 with a "green roof"). 

>

Int(Int(Int(f(x,y,z),z=0..1-1/3*y), y=3-3*x^2..3), x=0..1) + Int(Int(Int(f(x,y,z),z=0..x^2), y=0..3-3*x^2), x=0..1) = evalf(int(int(int(f(x,y,z),z=0..1-1/3*y), y=3-3*x^2..3), x=0..1) + int(int(int(f(x,y,z),z=0..x^2), y=0..3-3*x^2), x=0..1)); Int(Int(Int(f(x,y,z),z=0..1-1/3*y), x=sqrt(1-1/3*y)..1), y=0..3) + Int(Int(Int(f(x,y,z),z=0..x^2), x=0..sqrt(1-1/3*y)), y=0..3) = evalf(int(int(int(f(x,y,z),z=0..1-1/3*y), x=sqrt(1-1/3*y)..1), y=0..3) + int(int(int(f(x,y,z),z=0..x^2), x=0..sqrt(1-1/3*y)), y=0..3)); Int(Int(Int(f(x,y,z),y=0..3-3*z), z=0..x^2), x=0..1) = evalf(int(int(int(f(x,y,z),y=0..3-3*z), z=0..x^2), x=0..1)); Int(Int(Int(f(x,y,z),y=0..3-3*z), x=sqrt(z)..1), z=0..1) = evalf(int(int(int(f(x,y,z),y=0..3-3*z), x=sqrt(z)..1), z=0..1)); Int(Int(Int(f(x,y,z),x=sqrt(z)..1), z=0..1-1/3*y), y=0..3) = evalf(int(int(int(f(x,y,z),x=sqrt(z)..1), z=0..1-1/3*y), y=0..3)); Int(Int(Int(f(x,y,z),x=sqrt(z)..1), y=0..3-3*z), z=0..1) = evalf(int(int(int(f(x,y,z),x=sqrt(z)..1), y=0..3-3*z), z=0..1));

 

 

 

 

 

 

	(43)

 

Example: Let and consider the region which is above below , in front of , bounded on the left by , and on the right by .  

>	f := (x,y,z) -> x^3*y^2*z*exp(x+2*y+3*z): 'f(x,y,z)' = f(x,y,z);

 

(44)

 

    To plot the region of integration I will plot the surfaces which make up the sides of region.  

>	solve({y=x, y=2-x^2});

 

(45)

 

>

bottomPlot := plot3d(1,x=0..1,y=x..2-x^2,color=blue):   topPlot := plot3d(2,x=0..1,y=x..2-x^2,color=red):  backPlot := plot3d([0,y,z],y=0..2,z=1..2,color=gray):  leftPlot := plot3d([x,x,z],x=0..1,z=1..2,color=green): rightPlot := plot3d([x,2-x^2,z],x=0..1,z=1..2,color=yellow): display({bottomPlot,topPlot,backPlot,leftPlot,rightPlot},scaling=constrained,axes=boxed,viewpoint=circleleft);

 

Let's compute this integral in all 6 orders of integration.  

Note: dx dy dz and dx dz dy (thinking of as an simple region) will need to be split into 2 pieces (1 with from "gray" to "green" and 1 from "gray" to "yellow"). 

>	evalf(int(int(int(f(x,y,z),z=1..2), y=x..2-x^2), x=0..1));

 

(46)

 

>	evalf(int(int(int(f(x,y,z),z=1..2), x=0..y), y=0..1) + int(int(int(f(x,y,z),z=1..2), x=0..sqrt(2-y)), y=1..2));

 

(47)

 

>	evalf(int(int(int(f(x,y,z),y=x..2-x^2), z=1..2), x=0..1));

 

(48)

 

>	evalf(int(int(int(f(x,y,z),y=x..2-x^2), x=0..1), z=1..2));

 

(49)

 

>	evalf(int(int(int(f(x,y,z),x=0..y), z=1..2), y=0..1) + int(int(int(f(x,y,z),x=0..sqrt(2-y)), z=1..2), y=1..2));

 

(50)

 

>	evalf(int(int(int(f(x,y,z),x=0..y), y=0..1), z=1..2) + int(int(int(f(x,y,z),x=0..sqrt(2-y)), y=1..2), z=1..2));

 

(51)

 

Example: Consider the solid region E bounded by the surfaces and . Let's begin by graphing this region. Notice that the bottom of the region is an ellipsoid opening upward and the top is a horizontal plane. The intersection of these two surfaces is the ellipse. 

>	solve(x^2/4+y^2/9=1,y);

 

(52)

 

>	bottom := plot3d(x^2/4+y^2/9,x=-2..2,y=-3/2*sqrt(4-x^2)..3/2*sqrt(4-x^2)): top := plot3d(1,x=-2..2,y=-3/2*sqrt(4-x^2)..3/2*sqrt(4-x^2)): display({bottom,top},scaling=constrained,axes=boxed,orientation=[35,75,10]);

 

Let's compute the centriod of this region.  

>

m := int(int(int(1,z=x^2/4+y^2/9..1),y=-3*sqrt(1-x^2/4)..3*sqrt(1-x^2/4)),x=-2..2); Myz := int(int(int(x,z=x^2/4+y^2/9..1),y=-3*sqrt(1-x^2/4)..3*sqrt(1-x^2/4)),x=-2..2); Mxz := int(int(int(y,z=x^2/4+y^2/9..1),y=-3*sqrt(1-x^2/4)..3*sqrt(1-x^2/4)),x=-2..2); Mxy := int(int(int(z,z=x^2/4+y^2/9..1),y=-3*sqrt(1-x^2/4)..3*sqrt(1-x^2/4)),x=-2..2); Centroid := (Myz/m,Mxz/m,Mxy/m);

 

 

 

 

 

	(53)

 

Example:  Consider the region of integration where is bounded above by and below by .  

>

top := plot3d(sqrt(25-x^2-y^2), x=-3..3, y=-sqrt(9-x^2)..sqrt(9-x^2), color=red): bottom := plot3d(4, x=-3..3, y=-sqrt(9-x^2)..sqrt(9-x^2), color=green): sphereGhost := plot3d(5-0.05,theta=0..2*Pi,phi=0..Pi/2,coords=spherical,color=gray,style=line): display({bottom,top,sphereGhost},scaling=constrained,axes=boxed,orientation=[45,45,0]);

 

Let's compute .  First, we'll compute the integral in spherical coordinates (a natural choice). 

ρ ranges from the plane out to the sphere. The plane becomes or in spherical coordinates. The sphere becomes . Obviously θ ranges from 0 to 2π. Finally, φ ranges from 0 to the circle where the sphere and the plane intersect. So we have and . Thus and so or . Recall that . On the sphere we have so we get and thus . 

>	Int(Int(Int((rho^2)^(-1/2) * rho^2*sin(phi), rho=4*sec(phi)..5), phi=0..arcsin(3/5)), theta=0..2*Pi) = int(int(int((rho^2)^(-1/2) * rho^2*sin(phi), rho=4*sec(phi)..5), phi=0..arcsin(3/5)), theta=0..2*Pi);

 

(54)

 

Now let's compute this integral using cylidrical coordinates. First, ranges from the plane to the (upper-half of the) sphere . In the last part we established that the circle where the sphere and plane intersect is .  

>	Int(Int(Int((r^2+z^2)^(-1/2) * r, z=4..sqrt(25-r^2)), r=0..3), theta=0..2*Pi) = int(int(int((r^2+z^2)^(-1/2) * r, z=4..sqrt(25-r^2)), r=0..3), theta=0..2*Pi);

 

(55)

 

By the way, if you attempt to make Maple compute this integral directly in rectangular coordinates, Maple will choke. If you would like to see Maple choke and die, just remove the "#" (comment symbol) and hit enter. When working in rectangular coordinates, MAPLE GIVES THE WRONG ANSWER!!! 

>	Int(Int(Int((x^2+y^2+z^2)^(-1/2),z=4..sqrt(25-x^2-y^2)),y=-sqrt(9-x^2)..sqrt(9-x^2)),x=-3..3) = int(int(int((x^2+y^2+z^2)^(-1/2),z=4..sqrt(25-x^2-y^2)),y=-sqrt(9-x^2)..sqrt(9-x^2)),x=-3..3);

 

(56)

 

While we've gone to all this trouble setting up bounds and such, let's compute the centroid of our region of integration. 

>

vol := int(int(int(rho^2*sin(phi),rho=4*sec(phi)..5),phi=0..arccos(4/5)),theta=0..2*Pi); Myz := int(int(int(rho*cos(theta)*sin(phi)*rho^2*sin(phi),rho=4*sec(phi)..5),phi=0..arccos(4/5)),theta=0..2*Pi); Mxz := int(int(int(rho*sin(theta)*sin(phi)*rho^2*sin(phi),rho=4*sec(phi)..5),phi=0..arccos(4/5)),theta=0..2*Pi); Mxy := int(int(int(rho*cos(phi)*rho^2*sin(phi),rho=4*sec(phi)..5),phi=0..arccos(4/5)),theta=0..2*Pi); Centroid := [Myz/vol, Mxz/vol, Mxy/vol];

 

 

 

 

 

	(57)

 

Example: Consider the region of integration where is above the cone and inside the sphere .  

First, let's plot (the surface of) our region. It makes sense to represent the cone using cylindrical coordinates (note the appearance of ). We have . So parametrically the cone is The sphere is more easily described in spherical coordinates as , so parametrically we have . For both surfaces (obviously) . Also, 0 is a lower bound for both and . The upper bound for both of these variables is determined by where the cone and sphere intersect. So we need to solve the equations (the cone) and (the sphere) simultaneously.  

>	allvalues(solve({z=5*r,r^2+z^2=4}));

 

(58)

 

So the upper bound for is . To find the upper bound for φ we need to solve (the cone in spherical coordinates) and the sphere). 

>	solve({rho*cos(phi)=5*rho*sin(phi),rho=2});

 

(59)

 

Thus the upper bound for φ is . 

>	conePlot := plot3d([r*cos(theta),r*sin(theta),5*r],r=0..sqrt(2/13),theta=0..2*Pi, color=red): sphPlot := plot3d([2*cos(theta)*sin(phi),2*sin(theta)*sin(phi),2*cos(phi)],theta=0..2*Pi,phi=0..arctan(1/5), color=blue): display({sphPlot,conePlot},scaling=constrained,axes=boxed,orientation=[-40,60]);

 

Now let's compute using spherical coordinates. While fixing our plot above, we did all of the "hard work" finding the bounds. 

>	Int(Int(Int(rho^2*rho*cos(phi)*rho^2*sin(phi), rho=0..2), phi=0..arctan(1/5)), theta=0..2*Pi) = int(int(int(rho^2*rho*cos(phi)*rho^2*sin(phi), rho=0..2), phi=0..arctan(1/5)), theta=0..2*Pi);

 

(60)

 

Next, let's compute using cylidrical coordinates. Again, we already found the bounds while fixing our plot. 

>	Int(Int(Int((r^2+z^2)*z*r, z=5*r..sqrt(4-r^2)), r=0..sqrt(2/13)), theta=0..2*Pi) = int(int(int((r^2+z^2)*z*r, z=5*r..sqrt(4-r^2)), r=0..sqrt(2/13)), theta=0..2*Pi);

 

(61)

 

Example: Consider the region of integration where is above the -plane, inside the sphere and outside the cone . First we will plot the surfaces which bound the region . 

If we approach this problem using spherical and cylindrical coordinates, our plots and integrals will come out nice.  Let's where the sphere and cone intersect. To do this, we need to solve the equations: and simultaneously. In spherical coordinates we have (so and . 

>	solve(2*cos(phi)=2*sin(phi),phi);

 

(62)

 

In cylindrical coordinates, the intersection of the cone and sphere amounts to solving and . Thus and so . 

>

bottom := plot3d(0, x=-2..2, y=-sqrt(4-x^2)..sqrt(4-x^2), color=green): sphPlot := plot3d([2*cos(theta)*sin(phi),2*sin(theta)*sin(phi),2*cos(phi)],theta=0..2*Pi,phi=Pi/4..Pi/2, color=blue): conePlot := plot3d([r*cos(theta),r*sin(theta),r],r=0..sqrt(2),theta=0..2*Pi, color=red): display({bottom,sphPlot,conePlot},scaling=constrained,axes=boxed,orientation=[-40,60]);

 

Let's compute using spherical coordinates. All of the "hard work" was done while we were finding bounds for our plot. 

>	Int(Int(Int(rho^2*rho*cos(phi)*rho^2*sin(phi), rho=0..2), phi=Pi/4..Pi/2), theta=0..2*Pi) = int(int(int(rho^2*rho*cos(phi)*rho^2*sin(phi), rho=0..2), phi=Pi/4..Pi/2), theta=0..2*Pi);

 

(63)

 

Now let's compute using cylidrical coordinates. Again, we've already done the "hard work". 

Note: I recommend integrating with respect to r before z. If you want to do z first, you'll need to break this up into two pieces. 

>	Int(Int(Int((r^2+z^2)*z*r, r=z..sqrt(4-z^2)), z=0..sqrt(2)), theta=0..2*Pi) = int(int(int((r^2+z^2)*z*r, r=z..sqrt(4-z^2)), z=0..sqrt(2)), theta=0..2*Pi);

 

(64)

 

Example: Let's compute the centroid of the "bumpy sphere" which is described by the spherical coordinate equation where and . 

>	bumpyRho := 1 + (1/5)*sin(6*theta)*sin(5*phi); plot3d([bumpyRho*cos(theta)*sin(phi),bumpyRho*sin(theta)*sin(phi),bumpyRho*cos(phi)],theta=0..2*Pi,phi=0..Pi,numpoints=5000,viewpoint=circleleft);

 

Looking at the bumpy sphere we might guess that its centroid is the origin. Let's verify this by computing some integrals. Frist, we compute the mass (i.e. volume) by integrating over 1 (don't forget to include the Jacobian since we're working in spherical coordinates). Then we integrate over x for the moment about the yz-plane, then over y for the moment about the xz-plane, and finally over z for the moment about the xy-plane.  

>

m := int(int(int(1*rho^2*sin(phi),rho=0..bumpyRho),theta=0..2*Pi),phi=0..Pi); Myz := int(int(int((rho*cos(theta)*sin(phi))*rho^2*sin(phi),rho=0..bumpyRho),theta=0..2*Pi),phi=0..Pi); Mxz := int(int(int((rho*sin(theta)*sin(phi))*rho^2*sin(phi),rho=0..bumpyRho),theta=0..2*Pi),phi=0..Pi); Mxy := int(int(int((rho*cos(phi))           *rho^2*sin(phi),rho=0..bumpyRho),theta=0..2*Pi),phi=0..Pi);

 

 

 

 

	(65)

 

Then the centroid is: . So the centroid is origin (just as we guessed). 

What if we want the centroid of the upper-half of this solid? (Restrict φ to range from 0 to π/2.) 

>

top := plot3d([bumpyRho*cos(theta)*sin(phi),bumpyRho*sin(theta)*sin(phi),bumpyRho*cos(phi)],theta=0..2*Pi,phi=0..Pi/2,numpoints=2000): bottom := plot3d([rho*cos(theta)*sin(Pi/2),rho*sin(theta)*sin(Pi/2),rho*cos(Pi/2)],rho=0..subs(phi=Pi/2,bumpyRho),theta=0..2*Pi): display({top,bottom},scaling=constrained,numpoints=1000);

 

>

vol := int(int(int(rho^2*sin(phi), rho=0..bumpyRho), phi=0..Pi/2), theta=0..2*Pi); Myz := int(int(int(rho*cos(theta)*sin(phi) * rho^2*sin(phi), rho=0..bumpyRho), phi=0..Pi/2), theta=0..2*Pi); Mxz := int(int(int(rho*sin(theta)*sin(phi) * rho^2*sin(phi), rho=0..bumpyRho), phi=0..Pi/2), theta=0..2*Pi); Mxy := int(int(int(rho*cos(phi) * rho^2*sin(phi), rho=0..bumpyRho), phi=0..Pi/2), theta=0..2*Pi); Centroid := [Myz/vol, Mxz/vol, Mxy/vol];

 

 

 

 

 

	(66)

 

Example: Let be the trapezoidal region with vertices (1,0), (2,0), (0,1), and (0,2). Consider the following integral:. Let's compute this integral using the following change of variables:and . We can compute the Jacobian of this transform using the "Jacobian" command. 

>	InverseM, InverseJ := Jacobian([y-x,x+y],[x,y],determinant=true);

 

(67)

 

This isn't the correct Jacobian because we fed in the "new in terms of old" equations. We really, need to use equations which look like "x = ..." and "y = ...". Let's find these and compute the "right" Jacobian. 

>	solve({u=y-x,v=x+y},[x,y]);

 

(68)

 

>	M, J := Jacobian([-u/2+v/2,v/2+u/2],[u,v],determinant=true);

 

(69)

 

Of course in the end the Jacobian of the inverse transformation is just the inverse of the Jacobian! Oh well, I guess we really didn't have to compute it this way after all. 

I have plotted the region of integration (in the xy-plane) below. Transform this region using the change of variables and plot the corresponding region in the uv-plane. 

>	polygonplot([[1,0],[2,0],[0,2],[0,1]]);

 

Before plotting...notice we could calculate this integral directly, but it's a little messy since the bottom of the region is made up of 2 pieces: and . The top is given by: . 

>	Int(Int(cos((y-x)/(x+y)),y=1-x..2-x),x=0..1)+ Int(Int(cos((y-x)/(x+y)),y=0..2-x),x=1..2)= int(int(cos((y-x)/(x+y)),y=1-x..2-x),x=0..1)+ int(int(cos((y-x)/(x+y)),y=0..2-x),x=1..2);

 

(70)

 

>	evalf(int(int(cos((y-x)/(x+y)),y=1-x..2-x),x=0..1)+ int(int(cos((y-x)/(x+y)),y=0..2-x),x=1..2));

 

(71)

 

Back to plotting...Map each point (1,0), (2,0), etc using the equations and . I'll change the view so we can better see where the region lies. 

>	polygonplot([[-1,1],[-2,2],[2,2],[1,1]],view=[-2..2,0..2],color=green);

 

Now let's compute the integral using our change of variables. Notice the left hand side of the region is bounded by the line through (-2,2) and (-1,1) which is just . The right hand side is bounded by . Also, this region is more naturally thought of as a type II region since it's "bottom" is kind of complicated (it's made up of 3 pieces: , , and ). So let's integrate with respect to first and second. Finally, don't forget the Jacobian!! 

>	Int(Int(cos(u/v)*abs(-1/2),u=-v..v),v=1..2)= int(int(cos(u/v)*abs(-1/2),u=-v..v),v=1..2);

 

(72)

 

>	3/2*sin(1) = evalf(3/2*sin(1));

 

(73)

 

Therefore,  

[Notice that although we found the same answer without using the change of variables, it looked far more complicated and would have been a nightmare to compute "by hand".] 

Example: Let's make Maple compute the Jacobian for spherical coordinates. [The following command spits out a "Jacobian matrix" and the "Jacobian determinant". We just need the latter. Since the command spits out a sequence of answers, we put a sequence of labels on the left hand side of the assignment (i.e. "").] 

>	M,J := Jacobian([rho*cos(theta)*sin(phi),rho*sin(theta)*sin(phi),rho*cos(phi)],[rho,phi,theta],'determinant'); J := simplify(J);

 

 

	(74)

 

Example: In this example we will consider a 4D version of spherical coordinates. Let's call our new fourth spherical coordinate (the greek letter tau) and the fourth rectangular coordinate . 

,  ,  ,  and   where , , , and .   

>

x := (rho,theta,phi,tau) -> rho*cos(theta)*sin(phi)*sin(tau): 'x(rho,theta,phi,tau)' = x(rho,theta,phi,tau); y := (rho,theta,phi,tau) -> rho*sin(theta)*sin(phi)*sin(tau): 'y(rho,theta,phi,tau)' = y(rho,theta,phi,tau); z := (rho,theta,phi,tau) -> rho*cos(phi)*sin(tau): 'z(rho,theta,phi,tau)' = z(rho,theta,phi,tau); t := (rho,theta,phi,tau) -> rho*cos(tau): 't(rho,theta,phi,tau)' = t(rho,theta,phi,tau);

 

 

 

 

	(75)

 

Let's find the Jacobian matrix and determinant of this coordinate transformation. 

>	M,J := Jacobian([x(rho,theta,phi,tau),y(rho,theta,phi,tau),z(rho,theta,phi,tau),t(rho,theta,phi,tau)],[rho,theta,phi,tau],'determinant'); J := simplify(J);

 

 

	(76)

 

The 4D sphere (centered at the origin with radius ) is defined by . Let's rewrite this in spherical coordinates and then compute the hyper-volume enclosed by this sphere (i.e. the quadruple intergral of 1 over  ). 

>	sphere := x(rho,theta,phi,tau)^2+y(rho,theta,phi,tau)^2+z(rho,theta,phi,tau)^2+t(rho,theta,phi,tau)^2=R^2; sphere := simplify(sphere); Int(Int(Int(Int(1*J,rho=0..R),theta=0..2*Pi),phi=0..Pi),tau=0..Pi) = int(int(int(int(1*J,rho=0..R),theta=0..2*Pi),phi=0..Pi),tau=0..Pi);

 

 

 

	(77)

 

The area of a circle of radius is Differentiate the formula for the area with respect to and you get which is the circumference. The volume of a sphere of radius is Differentiate the formula for the volume with respect to and you get 4 which is the surface area! If we differentiate the "hypervolume" of your sphere, you will get the volume of its "surface". Check out Wikipedia: https://en.wikipedia.org/wiki/Sphere