An Extended Example of Triple Integration by Russell D. Blyth - blythrd@slu.edu and Michael K. May, S.J. - maymk@slu.edu restart; with(plottools): with(plots): This worksheet considers the following problem: Find the volume of the solid bounded by the planes z=-6, y=0, y-x=4 and 2x+y+z=4. The major challenge in this problem is to visualize the region so that appropriate limits of integration (and order of integration) can be determined.

We begin by plotting the four planes . We use the implicitplots3d command, color coding the 4 planes. eqn1 := z=-6; eqn2 := y=0; eqn3 := y-x=4; eqn4 := 2*x+y+z=4; planes := implicitplot3d([eqn1, eqn2, eqn3, eqn4], x=-5..6, y=-1..7, z=-7..15, color=[blue, red, green, yellow], style=patchnogrid, transparency=.5,axes=normal): display(planes); A first activity is to drag the plot around so that we look at it down each axis in turn. We want to note which axes give a nice front and back to the tetrahedron. If we look down the z-axis, the red and green planes become border lines, and the yellow and blue planes become front and back. If we look down the x-axis, the red and blue planes become border lines, and the yellow and green planes become front and back. If we look down the y-axis, the blue plane becomes a border line, the red plane becomes the back, and the yellow and green planes each become part of the front. This probably means we don't want to use the y-axis as the innermost integration. We also find the coordinates of the corners of the tetrahedron and add them to the picture. p1 := solve({eqn2, eqn3, eqn4},{x,y,z}); point1 := subs(p1,<x,y,z>); p2 := solve({eqn1, eqn3, eqn4},{x,y,z}): point2 := subs(p2,<x,y,z>); p3 := solve({eqn1, eqn2, eqn4},{x,y,z}): point3 := subs(p3,<x,y,z>); p4 := solve({eqn1, eqn2, eqn3},{x,y,z}): point4 := subs(p4,<x,y,z>); points := pointplot3d([point1, point2, point3, point4], symbol=diamond, symbolsize=30,color=gold,axes=normal): display([points,planes]); With the vertices known, we can create the edges of the tethedron as parameterized curves, using the fact that t*q1+(1-t)*q2 is the line segment from q1 to q2 as t goes from 0 to 1. line1 := t*point1 + (1-t)*point2: line2 := t*point1 + (1-t)*point3: line3 := t*point1 + (1-t)*point4: line4 := t*point2 + (1-t)*point3: line5 := t*point2 + (1-t)*point4: line6 := t*point3 + (1-t)*point4: linelist := map(convert, {line1, line2, line3, line4, line5, line6}, list): lines := spacecurve(linelist, t=0..1,color=black, thickness=2): display([points,planes, lines]); This shows the outline of the tetrahedron clearly.

To clean up the picture of the tetrahedron we replace the planes the triangular regions that make up exactly the faces of the tetrahedron. To do that we need to parameterize each face. We start with blue and red faces since they are easiest. For the blue face we look down the z-axis (orientation=[270,0]), whic shows the triangular section becomes a simple triangle in the plane. The three vertices on the plane are (2,6), (5,0), and (-4,0). eqn1; planeBlue := implicitplot3d([eqn1], x=-5..6, y=-1..7, z=-7..15, color=[blue], style=patchnogrid, transparency=.5,axes=normal): display([points, lines, planeBlue],axes=normal, orientation=[270,0]); The z value on the blue plane is always -6. The y values go from 0 to 6. For each y value the x values run from the green plane (eqn3) to the yellow plane (eqn4). We solve eqn3 and eqn4 for x in terms of y (with z=-6). solve(eqn3,x); solve(subs(z=-6,eqn4),x); Therefore, for each y, we see that x runs from y-4 to -y/2+5. We parameterize x for a given y with the same q1*t + q2*(1-t) method. s1 := plot3d([(y-4)*t+(-y/2+5)*(1-t),y,-6], t=0..1,y=0..6, transparency=.5, color=blue, style=patchnogrid): display([points, lines, s1],axes=normal); For the other three planes we look at the tetrahedron down the y axis (orientation=[270,90]). eqn2 := y=0; eqn3 := y-x=4; eqn4 := 2*x+y+z=4; otherplanes := implicitplot3d([eqn2, eqn3, eqn4], x=-5..6, y=-1..7, z=-7..15, color=[red, green, yellow], style=patchnogrid, transparency=.5,axes=normal): display([points, lines, otherplanes],axes=normal, orientation=[270,90]); In this view the red plane (eqn2) is the front face, the green face (eqn3) is the back face and the yellow face (eqn4) is the nonvertical plane. The lines of the tetrahedron form three triangles (the whole triangle, which is subdivided into a left and a right triangle). These triangles all have z running from -6 to 12. We now consider parametrizing each of these triangles one at a time. The outer triangle is on the front face (red plane with y=0). For each z, x runs from the intersection of the green face (eqn3) and the red face (eqn2) to the intersection of the yellow face (eqn4) and the red face (eq2). We solve the equations eqn3 and eqn4 for x in terms of z (with y=0). solve({subs(y=0,eqn3)},x); solve({subs(y=0,eqn4)},x); Therefore for each value of z, we have x going from x=-4 to x=-z/2+2. s2 := plot3d([-4*t+(-z/2+2)*(1-t), 0, z], t=0..1,z=-6..12, transparency=.5, color=red, style=patchnogrid): display([points, lines, s2],axes=normal); The left triangle is on the green face. For each z, x runs from the intersection of the green face (eqn3) and the red face (eqn2) to the intersection of the green face (eqn3) and the yellow face (eq4). We solve the quations eqn2 and eqn4 for x in terms of z (with y=x+4). solve({subs(y=x+4,eqn2)},x); solve({subs(y=x+4,eqn4)},x); Therefore for each value of z, we have x going from x=-4 to x=-3/z. s3 := plot3d([-4*t+(-z/3)*(1-t),-4*t+(-z/3)*(1-t)+4,z], t=0..1,z=-6..12,transparency=.5, color=green, style=patchnogrid): display([points, lines, s3],axes=normal); The right triangle is on the right face. For each z, x runs from the intersection of the yellow face (eqn4) and the green face (eqn3) to the intersection of the yellow face (eqn4) and the red face (eq2). We solve the quations eqn2 and eqn3 for x in terms of z (with y=-2x-z+4). solve({subs(y=-2*x-z+4,eqn3)},x); solve({subs(y=-2*x-z+4,eqn2)},x); Therefore for each value of z, we have x going from x=-z/3 to -z/2+2. s4 := plot3d([-z/3*t+(-z/2+2)*(1-t), -2*(-z/3*t+(-z/2+2)*(1-t)) -z+4, z], t=0..1,z=-6..12,transparency=.5, color=yellow, style=patchnogrid): display([points, lines, s4],axes=normal); We plot all four faces together: display([points, lines, s1, s2, s3, s4],axes=normal); We can also look at a slice at height z=z0 somewhere between -6 and 12. This slice is a combination of (z0+6)/18 parts the top vertex = (-4,0,12) and (12-z0)/18 parts the yellow triangle = <(y-4)*t+(-y/2+5)*(1-t),y,-6> (the face opposite the vertex). We look first at the slice at height z0=2. z0:=2; s5 := plot3d(<-4,0,12>*(z0+6)/18+<(y-4)*t+(-y/2+5)*(1-t),y,-6>*(12-z0)/18, t=0..1,y=0..6, transparency=0,color=black,grid=[20,20]): display([points, lines, s5],axes=normal); With a little work we can make this an animation of the slices. z0 := 'z0': slice := z0 -> plot3d(<-4,0,12>*(z0+6)/18+<(y-4)*t+(-y/2+5)*(1-t),y,-6>*(12-z0)/18, t=0..1,y=0..6, transparency=0,color=black,grid=[20,20],axes=normal): display([seq(display([points, lines, slice(ht), s1, s2, s3, s4]),ht=-6..12)], insequence=true,axes=normal); To view the animation, select the picture, reduce the frame rate to 2 per second, set the loop to repeat, and start the animation. It may be clearer to look at the animation without the sides. display([seq(display([points, lines, slice(ht)]),ht=-6..12)], insequence=true,axes=normal); For a fixed height z0 and value y0 of y, we can also show the range of integration for x (you may change z0 to anything from -6 to 12 and y0 to anything from 0 to 4-z0/3). To determine the range of x for a given z0 and y), note that x runs from the green plane (eqn3) to the yellow plane (eqn4) solve(eqn3,x); solve(eqn4,x); So x runs from y0-4 to 2 - y0/2 - z0/2. z0 := 2; y0:= 1.5; s5 := plot3d(<-4,0,12>*(z0+6)/18+<(y-4)*t+(-y/2+5)*(1-t),y,-6>*(12-z0)/18, t=0..1,y=0..6, transparency=0,color=black,grid=[20,20]): xr:=spacecurve([t*(y0-4)+(1-t)*(2-y0/2-z0/2),y0,z0],t=0..1,color=white,thickness=3): display([points, lines, s5,xr], axes=normal);

We can now set up the integral. We use the order of integration that matches the plots at the end of the previous section. For given values of y and z we integrate x from y-4 to (4-y-z)/2. I1 := Int(integrand,x=(y-4)..(4-y-z)/2); We look at the frame of our tetrahedron again, from the viewpoint of looking down the x axis. display(lines, orientation=[0,90], axes=normal, labels=["x", "y", "z"]); We see that we can next integrate in either y or z. If we integrate in y next (to match integrating over the black triangle above at any height z) we go from the line y=0 (the projection of the red plane on the y-z plane) to the line y=4-z/3 (the intersection of eq3 and eq4, the yellow and green planes, projected into the y-z plane). solve({eqn3,eqn4},{x,y}); I2 := Int(I1, y=0..(4-z/3)); That leaves us ready to integrate with respect to z. We see that z goes from -6 to 12. I3 := Int(I2,z=-6..12); To find the volume, we set integrand to 1. I31 := eval(I3,integrand=1); value(I31); We can also show this integral in steps: Int(Int(Int(1, x=y-4..2-y/2-z/2),y=0..4-z/3),z=-6..12) = Int(Int(int(1, x=y-4..2-y/2-z/2),y=0..4-z/3),z=-6..12); Int(Int(Int(1, x=y-4..2-y/2-z/2),y=0..4-z/3),z=-6..12) = Int(expand(int(int(1, x=y-4..2-y/2-z/2),y=0..4-z/3)),z=-6..12); Int(Int(Int(1, x=y-4..2-y/2-z/2),y=0..4-z/3),z=-6..12) = int(int(int(1, x=y-4..2-y/2-z/2),y=0..4-z/3),z=-6..12); To find the x (respectively, y and z) coordinate of the center of gravity of the tetrahedron we divide the integral with integrand x (respectively, y and z) by the integral with integrand 1. cgx := subs(integrand=x,I3)/eval(I3,integrand=1); value(cgx); cgy := subs(integrand=y,I3)/eval(I3,integrand=1); value(cgy); cgz := subs(integrand=z,I3)/eval(I3,integrand=1); value(cgz);