Contents

Evaluating Triple Iterated Integrals

In passing from double to triple integrals, there is much less that is novel than in passing from single to double integrals.
MATLAB has a built-in triple integrator triplequad similar to dblquad, but again, it only integrates over rectangular boxes. However, more general threefold iterated integrals can be formulated
and evaluated using the symbolic toolbox's routines, in the same manner as twofold iterated integrals. And if the first of
the three integrals can be evaluated symbolically, it is then possible to use dblquad for the remaining integrations.

Example 1

We will begin by considering the iterated integral

We can evaluate this easily as

syms xyz
firstans=int(int(int(x*y*z,z,0,3-x-y),y,0,3-x),x,0,3)

firstans =
81/80

We can also evaluate the first integral only, and then use dblquad or numint2 for the remaining two integrations:

f=int(x*y*z,z,0,3-x-y)
secondans=numint2(f,y,0,3-x,x,0,3)

f =
(x*y*(x + y - 3)^2)/2
secondans =
1.0125

Let's see how the answers compare:

double(firstans-secondans)

ans =
0

To use triplequad here, we would again have to make a change of variables

So in the new coordinates, the integral is over the box 0<x<3, 0<u<1, 0<v<1 and the integrand is

But because dblquad is sometimes slow (and the same algorithm extrapolated to three dimensions would be even slower) or doesn't work very well,
we have looked for an alternative, based on a freeware package of M-files called the numerical integration toolbox. This toolbox is much faster and more adaptable (though admittedly sometimes harder to use) than dblquad. Copy this package from http://www.math.umd.edu/users/jmr/241/mfiles/nit/ into a subdirectory of your working directory called nit. To use this package, first add the package with the command:

addpath nit

To see how it works, we first try out an alternative version of numint2, based on nit, called newnumint2. The syntax is the same as for numint2, but there is an optional final parameter, which is the number of subdivision points to use in each direction. (The default
is 20.) The higher this parameter, the more accurate the results; however, the running time of the algorithm is proportional
to the square of this parameter, so be careful about increasing it.

fourthans=newnumint2(f,y,0,3-x,x,0,3)

fourthans =
1.0125

Note that this answer looks the same as the first two to four decimal places, but that we are spared all the annoying warning
messages. The accuracy is:

double(firstans-fourthans)

ans =
0

which is pretty good!! Now we try a totally numerical integration. The analogue of newnumint2 for triple integrals is called numint3. Again there is an optional final parameter, which is the number of subdivision points to use in each direction. (The default
is 8.) The higher this parameter, the more accurate the results; however, the running time of the algorithm is proportional
to the cube of this parameter, so be extra careful about increasing it.

fourthans=numint3(x*y*z,z,0,3-x-y,y,0,3-x,x,0,3)

fourthans =
1.0125

The accuracy is:

double(firstans-fourthans)

ans =
0

Again pretty good!

Problem 1

Evaluate the threefold iterated integral

both symbolically and numerically. (For the numerical integration, first do the integral over z, then use newnumint2 on this intermediate result.)

Finding the Limits

The remaining issue in the evaluation of triple integrals is the determination of limits. Once this has been done, the actual
integration can always be performed by one of the methods illustrated above.

Example 2

Let us consider the region above the paraboloid

and inside the ellipsoid

We might begin by plotting the two surfaces simultaneously. Since z is never negative on or above the paraboloid, we are interested only in positive values of z. Consequently, we will solve the equation of the ellipsoid for positive z and use ezmesh to plot both surfaces.

From the viewpoint we have chosen, using MATLAB's view comand, we can see that there is a region inside the ellipsoid that is bounded below by the paraboloid and above by the ellipsoid.
This determines the lower and upper limits for z. Accordingly, we set

z1=2*x^2 + 3*y^2
z2=sqrt(20-3*x^2-2*y^2)

z1 =
2*x^2 + 3*y^2
z2 =
(20 - 2*y^2 - 3*x^2)^(1/2)

To determine the limits for x and y, it may be helpful to plot the curve along which the upper and lower limits for z coincide. We will do this using ezplot:

ezplot(z2-z1, [-3, 3, -2, 2])

Although this looks like an ellipse, it is actually a quartic curve. We will establish the limits by solving the equation
of the curve for y as a function of x. We expect two solutions, corresponding to the upper and lower y limits, and we should be able to determine the x limits by setting the y limits equal to each other.

Here solve has found four solutions instead of two. This should not surprise us too much since the equation in question is quartic.
Of the four solutions, however, two of them are complex, as we can see by substituting x=0.

double(subs(ylims,x,0))

ans =
1.1763
0.0000 + 1.2673i
-1.1763
-0.0000 - 1.2673i

Thus solution #3 should be the lower limit for y, solution #1 should be the upper limit, and solutions #2 and 4 are complex
and should be discarded. So we set

It takes a proposed set of limits of integration (written from the inside to the outside) and draws a picture of the corresponding
region.

viewSolid(z,z1,z2,y,y1,y2,x,x1,x2)

This looks reasonable; we may now use the limits we have found to integrate any function we like over the region in question.

Problem 2

Find the appropriate limits for integrating over the region bounded by the paraboloids

and

Check your results using viewSolid.

What is special about the examples we have just considered is that the region in question is bounded by two surfaces, each
of which has an equation specifying z as a function of x and y. The procedure we have used in this case will work for any such region, as long as solve is able to solve the associated
equations. We will look now at an example of another sort.

Example 3

We consider the region interior to both of the elliptical cylinders, x^2+2y^2 = 9 and 2x^2+z^2 = 8. We begin by observing that the condition that a point be interior to both cylinders constrains x^2 to be at most 4. We also observe that, for fixed x, the limits for y and z are determined independently of each other. So we'll want to put the x-integral on the outside. We can immediately set

x1=-2; x2=2;
y1=-sqrt((9-x^2)/2); y2=-y1;
z1=-sqrt(8-2*x^2); z2=-z1;

We can again use viewSolid to check the plausibility of these limits, although the resulting plot is more difficult to interpret in this case, since
it shows only the upper and lower boundaries as surfaces; the side boundaries are shown as ensembles of vertical lines.