VC.08 Transforming 3D Integrals
 BASICS

\(Off[General::spell]; 
Off[General::spell1]; 

    \(Off[ParametricPlot3D::ppcom]; \)\)Off[ParametricPlot::ppcom]; Off[Plot::plnr]; SetOptions[Limit, Analytic -> True];

CMView = {2.7, 1.6, 1.2};

Vector Drawer

BeginPackage["Vector3D`", "Graphics`Colors`"];

\(Axes3D::"usage" = 
      "Axes3D[a, b] creates a Graphics3D object of cartesian axes with x, y, and z running from -a/3 to a, and with axes labels b units beyond the tips of the axes.  Axes3D[a] is Axes3D[a, a/8]."; 

    \(Perpend::"usage" = 
      "Perpend[a] returns a unit vector which is perpendicular to the vector a."; \)\)Arrow::"usage" = 
      "Arrow[a] produces a 2 or 3 dimensional vector from the origin to a. \
Arrow[a, Tail -> tail] gives a vector from tail to a."; Vector::"usage" = 
      "Vector[a, b] gives a 2D or 3D vector from a to b."; ArrowHead::"usage" = 
      "ArrowHead[a, vec] produces an arrowhead with its tip placed at the \
point a, pointing in the direction of the vector vec."; Tail::"usage" = 
      "Tail -> point puts the tail of the vector at point."; Norm::"usage" = 
      "Norm[vec] returns the Euclidean length of vec."; Normalize::"usage" = 
      "Normalize[vec] returns a unit vector in the direction of vec, \
provided vec is nonzero."; Aperture::"usage" = 
      "Aperture is the ratio of the radius of the base to the length of the \
head of a vector."; TipSize::"usage" = 
      "TipSize specifies an absolute size for the inner tip length of the \
head of a vector."; TipRatio::"usage" = 
      "TipRatio is the ratio of the inner length of the arrowhead to the \
length of the shaft of a vector."; HeadRatio::"usage" = 
      "HeadRatio is the ratio of the outer length of the arrowhead to the \
length of the shaft of a vector."; HeadSize::"usage" = 
      "HeadSize specifies an absolute size for the head of a vector."; 
    ScaleFactor::"usage" = 
      "ScaleFactor specifies the amount to scale a vector in length and may \
be either a number or a pure function."; ZeroVectorPointSize::"usage" = 
      "ZeroVectorPointSize is the size of the point used to represent a \
zero vector."; Options[Vector] = 
      Options[Arrow] = 
        Options[ArrowHead] = {HeadRatio -> 0.18, 
            TipRatio -> 0.14, Aperture -> 0.3, 
            ScaleFactor -> 1, HeadType -> Polygon, 
            ShaftQ -> True, EdgesQ -> True, 
            VectorColor -> RGBColor[0, 0, 1], ShaftWidth -> 0.005, 
            ZeroVectorPointSize -> 0.01}; Begin["`Private`"]; SetOptions[ParametricPlot, AspectRatio -> Automatic]; SetOptions[Plot, AspectRatio -> Automatic]; SetOptions[Graphics, AspectRatio -> Automatic]; Axes3D[u_, v_] := 
      Graphics3D[{{Blue, Line[{{-(u/3), 0, 0}, {u, 0, 0}}]}, 
          Text["x", {u + v, 0, 0}], {Blue, 
            Line[{{0, -(u/3), 0}, {0, u, 0}}]}, 
          Text["y", {0, u + v, 0}], {Blue, 
            Line[{{0, 0, -(u/3)}, {0, 0, u}}]}, 
          Text["z", {0, 0, u + v}]}]; Axes3D[u_] := Axes3D[u, (u/8)]; Norm[a_?VectorQ] := N[(Sqrt[a . a])]; Normalize[a_?VectorQ] := If[a . a == 0, a, (a/Norm)[a]]; Perpend[a_?VectorQ] := 
      Normalize[{a[[2]], 
            -a[[1]]}] /; 
        Length[a] == 2Perpend[a_?VectorQ] := 
      Normalize[
          If[a[[1]] == 0, {1, 0, 0}, {0, 
              a[[3]], 
              -a[[2]]}]] /; 
        Length[a] == 3Base3D = 
      Table[{Cos[((2 k π)/8.) ], Sin[((2 k π)/8.) ]}, {k, 0, 
          8}]; Vector[from_?VectorQ, to_?VectorQ, opts___] := 
      Module[{sq, sw, ht, hs, hr, tr, ap, sf, vc, zvps, diff = to - from, 
            tip}, {sq, sw, ht, hs, hr, ts, tr, ap, sf, vc, zvps} = 
            {ShaftQ, ShaftWidth, HeadType, HeadSize, HeadRatio, TipSize, 
                  TipRatio, Aperture, ScaleFactor, VectorColor, 
                  ZeroVectorPointSize} /. {opts} /. Options[Vector]; 
          If[from == to, 
            Return[Graphics[{PointSize[zvps], vc, Point[from]}]]]; 
          diff = If[NumberQ[N[sf]], sf diff, sf[diff]]; tip = from + diff; 
          If[NumberQ[hs], hr = (hs/Norm)[diff]; tr = 0.75 hr]; 
          If[NumberQ[ts], tr = (ts/Norm)[diff]]; 
          Graphics[{
              If[sq, {vc, Thickness[sw], Line[{from, tip - tr diff}]}, {}], 
              With[{trans = ap hr Norm[diff] Perpend[diff]}, {vc, 
                  Thickness[sw], 
                  ht[{tip, tip - hr diff + trans, tip - tr diff, 
                      tip - hr diff - trans, tip}]}]}]] /; 
        Length[from] == Length[to] == 2Vector[from_?VectorQ, to_?VectorQ, opts___] := 
      Module[{sq, eq, sw, hs, hr, tr, ap, sf, vc, zvps, diff = to - from, 
            tip}, {sq, eq, sw, hs, hr, ap, sf, vc, zvps} = 
            {ShaftQ, EdgesQ, ShaftWidth, HeadSize, HeadRatio, Aperture, 
                  ScaleFactor, VectorColor, ZeroVectorPointSize} /. {opts}
                   /. Options[Vector]; 
          If[from == to, 
            Return[Graphics3D[{PointSize[zvps], vc, Point[from]}]]]; 
          diff = If[NumberQ[N[sf]], sf diff, sf[diff]]; tip = from + diff; 
          If[NumberQ[hs], hr = (hs/Norm)[diff]]; tr = (hr/2); 
          Graphics3D[{
              If[sq, {vc, Thickness[sw], Line[{from, tip - tr diff}]}, {}], 
              If[eq, {}, EdgeForm[]], {SurfaceColor[vc], 
                With[{temp = Perpend[diff]}, 
                  (Polygon[Append[#1, tip]]&)/@
                    Partition[
                      (tip - hr diff + #1&)/@
                        ((ap hr Norm[diff] Base3D) . {temp, 
                              Normalize[temp\[Cross]diff]}), 2, 1]]}}]] /; 
        Length[from] == Length[to] == 3Arrow[a_, Tail -> b_, opts___] := Vector[b, a + b, opts]Arrow[a_, opts___] := Vector[Table[0, {Length[a]}], a, opts]\(ArrowHead[a_, b_, opts___] := 
      Vector[a - b, a, ShaftQ -> False, opts]; 

End[]; 

    \(EndPackage[]; \)\)

Graphics`Colors`GosiaGreen = RGBColor[0,  0.392187,  0];

Spherical Setup

\(Clear[polar]; 

    polar[{x_, y_, z_}] := {\@\(x\^2 + y\^2 + z\^2\), ArcTan[y\/x], 
        ArcCos[z\/\@\(x\^2 + y\^2 + z\^2\)]}; 
Clear[X, length]; 

    length[X_] := \@\(X . X\); 
\(rr = 3; \)\)R = 5; scale = 0.8; red = RGBColor[1, 0, 0]; blue = RGBColor[0, 0, 1]; color1 = RGBColor[0.8, 0.2, 0.8]; color2 = RGBColor[0.2, 0.8, 0.2]; black = RGBColor[0, 0, 0]; \(Clear[axes]; 

    axes := Graphics3D[{Line[{X, {1.5\ rr, 0, 0}}], 
          Line[{X, {0, 1.5\ rr, 0}}], Line[{X, {0, 0, 2\ rr}}]}]; 

    Clear[xyplane]; 

    xyplane := 
      Graphics3D[
        Polygon[{{\(-2\)\ rr, \(-2\)\ rr, 0}, {\(-2\)\ rr, 2\ rr, 0}, {2\ rr, 
              2\ rr, 0}, {2\ rr, \(-2\)\ rr, 0}}]]; 
\(X = {0, 0, 0}; \)\)origin = Graphics3D[{black, Point[X]}]; Z = {2, 3, 5}; \(Clear[point]; 

    point := Graphics3D[{red, PointSize[0.03], Point[scale\ Z]}]; 

    \(Z0 = {0, 0, length[0.8\ Z]}; \)\)\(Clear[zpoint]; 
zpoint = Graphics3D[{black, Point[Z0]}]; 

    \(XY = scale\ Z\ {1, 1, 0}; \)\)\(Clear[xypoint]; 
xypoint := Graphics3D[{black, Point[XY]}]; 

    \(X0 = {length[scale\ Z\ {1, 1, 0}], 0, 0}; \)\)xpoint = Graphics3D[{black, Point[X0]}]; Y0 = {0, length[scale Z {1, 1, 0}], 0}; ypoint = Graphics3D[{black, Point[Y0]}]; \(Clear[ray]; 
ray := Graphics3D[{blue, Line[{X, scale\ Z}]}]; 

    Clear[ray1]; 
ray1 := Graphics3D[{color1, Line[{scale\ Z, Z}]}]; 

    Clear[xyray]; 
xyray := Graphics3D[{red, Line[{X, Z\ {1, 1, 0}}]}]; 

    \(rays := {ray, xyray}; \)\)xyaux = 
      (length[0.8 Z] (0.8 Z {1, 1, 0}))/
        length[0.8 Z {1, 1, 0}]; \(Clear[zarc]; 

    zarc = Graphics3D[{blue, 
          Line[Table[
              \(Cos[t]\ \) {0, 0, length[0.8\ Z]} + Sin[t]\ xyaux, {t, 0, 
                \(polar[0.8\ Z]\)[[3
                  ]],  .1}]]}]; 
Clear[zarc1]; 

    zarc1 = Graphics3D[{color1, 
          Line[Table[
              \(Cos[t]\ \) {0, 0, length[0.8\ Z]} + Sin[t]\ xyaux, {t, 
                \(polar[0.8\ Z]\)[[3]], 
                1,  .1}]]}]; 
Clear[zarc2]; 

    zarc2 = Graphics3D[{color1, 
          Line[Table[
              \(Cos[t]\ \) {0, 0, length[0.8\ Z]} + Sin[t]\ xyaux, {t, 
                \(-0.8\), 0,  .1}]]}]; 
Clear[xyarc]; 

    xyarc = Graphics3D[{blue, 
          Line[Table[
              Cos[t]\ X0 + Sin[t]\ Y0, {t, 0, 
                \(polar[XY]\)[[2
                  ]],  .1}]]}]; 
Clear[xyarc1]; 

    xyarc1 = Graphics3D[{color1, 
          Line[Table[Cos[t]\ X0 + Sin[t]\ Y0, {t, \(- .5\), 0,  .1}]]}]; 

    Clear[xyarc2]; 

    xyarc2 = Graphics3D[{color1, 
          Line[Table[
              Cos[t]\ X0 + Sin[t]\ Y0, {t, 
                \(polar[XY]\)[[2]], 
                \(polar[XY]\)[[2]] + 
                  0.3,  .1}]]}]; 
Clear[stick]; 

    stick = Graphics3D[{red, Line[{scale\ Z, XY}]}]; 

    \(xaxislabel = Text["x", {0.363, 0.26}]; \)\)yaxislabel = Text["y", {0.73, 0.3}]; zaxislabel = Text["z", { .5, 0.67}]; stickplot = 
      Show[xyplane, axes, ray, point, Boxed -> False, 
        ViewPoint -> {2.747, 1.594, 1.168}, 
        Epilog -> {xaxislabel, yaxislabel, zaxislabel}, 
        DisplayFunction -> Identity]; rho = Text["r", {0.53, 0.45}]; rlabel = 
      Show[xyplane, axes, ray, point, Boxed -> False, 
        ViewPoint -> {2.747, 1.594, 1.168}, 
        Epilog -> {rho, xaxislabel, yaxislabel, zaxislabel}, 
        DisplayFunction -> Identity]; phi = Text["s", {0.55, 0.57}]; slabel = 
      Show[xyplane, zpoint, axes, ray, point, zarc, zarc1, zarc2, 
        Boxed -> False, ViewPoint -> {2.747, 1.594, 1.168}, 
        Epilog -> {phi, xaxislabel, yaxislabel, zaxislabel}, 
        DisplayFunction -> Identity]; theta = Text["t", { .5, 0.3}]; tlabel = 
      Show[xyplane, axes, xypoint, xpoint, ray, stick, point, xyarc, xyarc1, 
        xyarc2, Boxed -> False, ViewPoint -> {2.747, 1.594, 1.168}, 
        Epilog -> {theta, xaxislabel, yaxislabel, zaxislabel}, 
        DisplayFunction -> Identity]; all = 
      Show[xyplane, axes, ray, stick, point, zpoint, xpoint, xypoint, xyarc, 
        xyarc1, xyarc2, zarc, zarc1, zarc2, Boxed -> False, 
        ViewPoint -> {2.747, 1.594, 1.168}, 
        Epilog -> {theta, phi, rho, xaxislabel, yaxislabel, zaxislabel}, 
        DisplayFunction -> Identity];

Crisscross

\(Clear[x]; 

    crisscross = 
      Plot[{Cos[x], Sin[x]}, {x, π\/8, \(3\ π\)\/8}, 
        PlotStyle -> {{RGBColor[1, 0, 0], Thickness[0.02]}}, 
        Axes -> None, AspectRatio -> 1\/3, 
        DisplayFunction -> Identity]; \)

B.1)  3D integrals

B.1.a)

Here's a plot of the solid R whose top skin is the surface
       z=4- (x^2)- ((y^2)/4)
and whose bottom skin lies on the plane z = 1.

topskin = 
      Plot3D[4 - (x^2) - ((y^2)/4), {x, -1.8, 1.8}, {y, -3.5, 3.5}, 
        PlotRange -> {0.9, 5}, ClipFill -> None, 
        DisplayFunction -> Identity]; \(base = 
      Graphics3D[
        Polygon[{{\(-2\), \(-4\), 1}, {2, \(-4\), 1}, {2, 4, 1}, {\(-2\), 4, 
              1}}]]; 


    solid = Show[topskin, base, Boxed -> False, 
        BoxRatios -> Automatic, AxesLabel -> {"x", "y", "z"}, 
        ViewPoint -> CMView, DisplayFunction -> $DisplayFunction]; 
    \)

You could use a 2D integral to measure the volume of this solid. But just to get experience, use a 3D integral to measure the volume of this solid.

Just as 
       (∫_^)(∫_( R)^) 
          dx dy
measures the area of a region R in two dimensions,
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
            dx dy dz\)\)\)
measures the volume of a solid R in three dimensions.
You are free to choose any order you like for the variables so that
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
              dx dy dz\)\) = 
        \(∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
                dz dy dx\)\) = 
          ∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
                dz dx 
                  dy\)\)\)\).
In the current situation, the top is
       z=4- (x^2)- ((y^2)/4)
and the bottom is
       z = 1
With z so conveniently displayed, it's natural to integrate with respect to z first. 
To set the limits for the first integral, fix x and y and enter the lowest and the highest values of z for this fixed x and y: 
       (∫_(zlow[x, y])^(zhigh[x, y])) dz.

\(Clear[zhigh, zlow, x, y, z]; 
    
{zlow[x_, y_], zhigh[x_, y_]} = {1, 4 - x\^2 - y\^2\/4}; 

    firstintegral = 
      ∫\_\(zlow[x, y]\)\%\(zhigh[x, y]\)\[DifferentialD]z\)

For the second integral, you've got your choice of whether to integrate with respect to y or x next.  Go with x.
The main problem is deciding what limits to insert for xhigh[y] and xlow[y]  in the second integral
       (∫_(xlow[y])^(xhigh[y]
            ))(∫_(zlow[x, y])^(zhigh[x, y])) 
          dz dx.
The key to deciding how to set xhigh[y] and xlow[y] is to ask:
For a fixed y, what are the lowest and the highest values x can have no matter what z is?  
The simplest way of answering this is to plot the shadow of this solid on xy-paper:

shadoweqn = Eliminate[{z == zhigh[x, y], z == zlow[x, y]}, z]

shadowcurves = Solve[shadoweqn, x]

\(Clear[xlow, xhigh]; 
xlow[y_] = \(-\(\@\(12 - y\^2\)\/2\)\); 

    xhigh[y_] = \@\(12 - y\^2\)\/2; 

    shadowplot = 
      Plot[{xlow[y], xhigh[y]}, {y, \(-2\) \@3, 2 \@ 3}, 
        PlotStyle -> {Blue, Red}, AxesLabel -> {"y", "x"}]; 
    \)

The second integral (∫_(xlow[y])^(xhigh[y]
            ))(∫_(zlow[x, y])^(zhigh[x, y])) 
          dz dx is:

secondintegral = 
      (∫_(xlow[y])^(xhigh[y]))firstintegral 
          \[DifferentialD]x

For the third integral, you have no choice but to integrate with respect to y.  The main problem is deciding what limits to insert for yhigh and ylow in the third integral
       (∫_ylow^yhigh)
          (∫_(xlow[y])^(xhigh[y]
              ))(∫_(zlow[x, y])^(zhigh[x, y])) 
            dz dx dy.

The key to deciding how to set ylow and yhigh is to ask:
What are the lowest and the highest values y can have no matter what z and x are?  
The simplest way of getting a hold on ylow and yhigh is to look at the shadow plot again:

Show[shadowplot];

In this problem, ylow and yhigh happen where xlow[y] and xhigh[y] meet:

ends = Solve[xhigh[y] == xlow[y]]

{ylow, yhigh} = {-2 (Sqrt[3]), 2 (Sqrt[ 3])}

The volume is measured by the third integral
       (∫_ylow^yhigh)
          (∫_(xlow[y])^(xhigh[y]
              ))(∫_(zlow[x, y])^(zhigh[x, y])) 
            dz dx dy:

volume = (∫_ylow^yhigh) secondintegral \[DifferentialD]y

Or:

NIntegrate[secondintegral, {y, ylow, yhigh}]

You can get this with a single 3D integral instruction for the volume measurement
       (∫_ylow^yhigh)
          (∫_(xlow[y])^(xhigh[y]
              ))(∫_(zlow[x, y])^(zhigh[x, y])) 
            dz dx dy:

threeDintegral = 
      (∫_ylow^yhigh)
          (∫_(xlow[y])^(xhigh[y]
              ))(∫_(zlow[x, y])^(zhigh[x, y]))1 
              \[DifferentialD]z \[DifferentialD]x \[DifferentialD]y

Good.
However, doing this in one step does not relieve you of the experience of setting zhigh[x, y], zlow[x, y], xhigh[y], xlow[y], yhigh, and ylow.
Tough break.

B.1.a.ii)

Discuss the physical meaning of the three integrals that were calculated in part i):

       first integral = 
        (∫_(zlow[x, y])^(zhigh[x, y]))dz, 
     
       second integral = 
        (∫_(xlow[y])^(xhigh[y]
              ))(∫_(zlow[x, y])^(zhigh[x, y])) 
            dz dx, and
     
       third integral = 
        (∫_ylow^yhigh)
            (∫_(xlow[y])^(xhigh[y]
                ))(∫_(zlow[x, y])^(zhigh[x, y])) 
              dz dx dy.

The first integral,
       (∫_(zlow[x, y])^(zhigh[x, y]))dz,
fixes x and y and measures the length of the stick that runs from {x, y, zlow[x, y]} to {x,  y,  zhigh[x, y]}. 
For example:

{x0, y0} = {0.5, 1}; stickinslice = 
      Graphics3D[{Blue, 
          Line[{{x0, y0, zlow[x0, y0]}, {x0, y0, zhigh[x0, y0]}}]}]; \(threedims = Axes3D[4.5]; 


    Show[stickinslice, base, threedims, BoxRatios -> Automatic, 
      Boxed -> False, ViewPoint -> CMView]; \)

The second integral
       (∫_(xlow[y])^(xhigh[y]
            ))(∫_(zlow[x, y])^(zhigh[x, y])) 
          dz dx
holds y = (Subscript[y, 0]) fixed but releases x and measures the accumulated area of the slice of the solid swept out by the sticks as x runs from xlow[y] to xhigh[y]:

\(sticksinslice = 
      Table[Graphics3D[{Blue, 
            Line[{{x, y0, zlow[x, y0]}, {x, y0, zhigh[x, y0]}}]}], {x, 
          xlow[y0], xhigh[y0], 1\/50\ \((xhigh[y0] - xlow[y0])\)}]; 

    \(Show[sticksinslice, base, threedims, Boxed -> False, 
      ViewPoint -> CMView]; \)\)

The third integral 
       (∫_ylow^yhigh)
          (∫_(xlow[y])^(xhigh[y]
              ))(∫_(zlow[x, y])^(zhigh[x, y])) 
            dz dx dy
releases y and measures the accumulated volume swept out by the slices as y runs from ylow to yhigh:

\(constituentslices = 
      Table[Graphics3D[{Blue, 
            Polygon[Table[{x, y, zhigh[x, y]}, {x, xlow[y] - 0.01, 
                  xhigh[y] + 0.01, 
                  1\/20\ \((xhigh[y] - xlow[y] + 0.02)\)}]]}], {y, ylow, 
          yhigh, \(yhigh - ylow\)\/20}]; 

    \(accumulatedslices = 
      Show[constituentslices, base, threedims, Boxed -> False, 
        ViewPoint -> CMView]; \)\)

The third integral measures the total volume.

B.1.a.iii) Triple integrals

Calculate 
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
            z dx dy dz\)\)\)
where R is the solid plotted in part i)

The layout is essentially the same as the calculation of the volume of R done in part a.i) above.

\(Clear[f]; 
f[x_, y_, z_] = z; \)

firstintegral = 
      (∫_(zlow[x, y])^(zhigh[x, y]))f[x, y, z] 
          \[DifferentialD]z

secondintegral = 
      (∫_(xlow[y])^(xhigh[y]))firstintegral 
          \[DifferentialD]x

thirdintegral = 
      (∫_ylow^yhigh) secondintegral \[DifferentialD]y

You can get this with one instruction for calculating
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
            f[x, y, z] dx dy 
              dz\)\)\)
        = 
        (∫_ylow^yhigh)
            (∫_(xlow[y])^(xhigh[y]
                ))(∫_(zlow[x, y])^(zhigh[x, y])) 
              f[x, y, z] dz dx 
                dy:

onestep = 
      (∫_ylow^yhigh)
          (∫_(xlow[y])^(xhigh[y]
              ))(∫_(zlow[x, y])^(zhigh[x, y]))f[x, y, z] 
              \[DifferentialD]z \[DifferentialD]x \[DifferentialD]y

Good.
Doing this in one step does not relieve you of the burden of setting 
       zhigh[x, y], zlow[x, y], 
       xhigh[y], xlow[y], 
       yhigh, and ylow.
You set them up for calculating
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
            f[x, y, z] dx dy 
              dz\)\)\)
the same way you set them up for calculating the volume of R by calculating
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
            dx dy dz\)\)\)
in part i).

B.2)  Transforming 3D integrals

B.2.a)

In two dimensions, you know that
       (∫_^)(∫_( (Subscript[R, xy]))^) 
          f[x, y] dx dy
        = 
        (∫_^)(∫_( (Subscript[R, uv]))^) 
            f[x[u, v], y[u, v]] (Subscript[A, xy])[u, v] du 
              dv.
where the area conversion factor (Subscript[A, xy])[u, v] is given by:

\(Clear[x, y, u, v, gradx, grady, Axy]; 

    gradx[u_, v_] = {D[x[u, v], u], D[x[u, v], v]}; 

    grady[u_, v_] = {D[y[u, v], u], D[y[u, v], v]}; 

    Axy[u_, v_] = Abs[Det[{gradx[u, v], grady[u, v]}]]\)

In three dimensions, what is the volume conversion factor
       (Subscript[V, xyz])[u, v, w] 
that makes
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\_xyz\)\%\ 
            f[x, y, z] dx dy 
              dz\)\)\)
       \(\( = 
        ∫\_\%\(∫\_\%\(∫\_\(\ R\_uvw\)\%\ 
              f[x[u, v, w], y[u, v, w], z[u, v, w]] V\_xyz[u, v, w] 
                du dv 
                dw\)\)\)\)

In mathematics, there is nothing like the power of analogy, and analogy works beautifully here. 
In two dimensions the area conversion factor (Subscript[A, xy])[u, v] is given by:

\(Clear[x, y, u, v, gradx, grady, Axy]; 

    gradx[u_, v_] = {D[x[u, v], u], D[x[u, v], v]}; 

    grady[u_, v_] = {D[y[u, v], u], D[y[u, v], v]}; 

    Axy[u_, v_] = Abs[Det[{gradx[u, v], grady[u, v]}]]\)

Analogously, in three dimensions the volume conversion factor 
(Subscript[V, xyz])[u, v, w] is given by:

\(Clear[x, y, z, u, v, gradx, grady, gradz, Vxyz]; 

    gradx[u_, v_, w_] = {D[x[u, v, w], u], D[x[u, v, w], v], 
        D[x[u, v, w], w]}; 

    grady[u_, v_, w_] = {D[y[u, v, w], u], D[y[u, v, w], v], 
        D[y[u, v, w], w]}; 

    gradz[u_, v_, w_] = {D[z[u, v, w], u], D[z[u, v, w], v], 
        D[z[u, v, w], w]}; 

    Vxyz[u_, v_, w_] = 
      Abs[Det[{gradx[u, v, w], grady[u, v, w], gradz[u, v, w]}]]\)

You just arrange the three gradients as (horizontal) rows of a matrix which most folks call the Jacobian matrix.

gradmatrix = {gradx[u, v, w], grady[u, v, w], gradz[u, v, w]}

Look at it in matrix form:

MatrixForm[gradmatrix]

Take the absolute value of its determinant:

Abs[Det[gradmatrix]]

This is the same as the volume conversion factor 
       (Subscript[V, xyz])[u, v, w]:

Vxyz[u, v, w]

Quite a pill to swallow, but who cares?  The computer is the one that has to swallow it.

B.2.b.i)

Check out the formula for the volume conversion factor 
       (Subscript[V, xyz])[u, v, w]
        = 
        |Det[{gradx[u, v, w], grady[u, v, w], 
              gradz[u, v, w]}]|
by using it to relate the volume of the ellipsoid
       ((((x/a)))^2)+((((y/b)))^2)+((((z/c)))^2)=(r^2) 
to the known volume ((4 π (r^3))/3) of the sphere of radius r
     (u^2) +  (v^2) +  (w^2) = (r^2).

\(Clear[x, y, z, a, b, c, r]; 

    ellipsoideqn = \((x\/a)\)\^2 + \((y\/b)\)\^2 + \((z\/c)\)\^2 == r\^2\)

Go to uvw-space with:

\(Clear[u, v, w]; 
    
{x[u_, v_, w_], y[u_, v_, w_], z[u_, v_, w_]} = {a\ u, b\ v, c\ w}; 
    \)

Here is the uvw-equation of the original ellipsoid:

ellipsoideqn /. {x -> x[u, v, w], y -> y[u, v, w], 
        z -> z[u, v, w]}

This tells you that if (Subscript[R, xyz]) is everything inside and on the ellipsoid
       ((((x/a)))^2)+((((y/b)))^2)+((((z/c)))^2)=(r^2) ,
then (Subscript[R, uvw]) is everything inside and on the sphere
     (u^2) +  (v^2) +  (w^2) = (r^2).
The volume of the sphere (Subscript[R, uvw]) is ((4 π (r^3))/3).
The  uvw-to-xyz-volume conversion factor is:

\(Clear[u, v, gradx, grady, gradz, Vxyz]; 

    gradx[u_, v_, w_] = {D[x[u, v, w], u], D[x[u, v, w], v], 
        D[x[u, v, w], w]}; 

    grady[u_, v_, w_] = {D[y[u, v, w], u], D[y[u, v, w], v], 
        D[y[u, v, w], w]}; 

    gradz[u_, v_, w_] = {D[z[u, v, w], u], D[z[u, v, w], v], 
        D[z[u, v, w], w]}; 

    Vxyz[u_, v_, w_] = 
      Abs[Det[{gradx[u, v, w], grady[u, v, w], gradz[u, v, w]}]]\)

So the xyz-volume measurement of (Subscript[R, xyz]) is
       a b c times the uvw-volume measurement of (Subscript[R, uvw]).
In short, the volume measurement of the ellipsoid
       ((((x/a)))^2)+((((y/b)))^2)+((((z/c)))^2)=(r^2) 
is 
       a b c((4  π  (r^3))/3).
Makes sense when you think about it.

B.2.b.ii)

Now get a bit more serious. 
Calculate 
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\_xyz\)\%\ 
            \(e\^\(-x\)\) Sin[y] dx dy 
              dz\)\)\)
where (Subscript[R, xyz]) is the region inside and on the box bounded by the planes
       z = 2  x -  1  (bottom).
       z = 2  x +  3  (top),
       y = x       (side),
       y = x +  4   (side),
       y = -2 x     (side), and
       y = -2 x +  5  (side).

A direct calculation would be a revolting task.  Even the task of plotting (Subscript[R, xyz]) wouldn't be pleasant.  But after you transform the integral to a favorable uvw-space, then everything is fairly easy.
Look at the formulas for the boundary surfaces:
       z = 2  x -  1 (bottom),  z = 2  x +  3 (top),
       y = x (side),  y = x +  4 (side),
       y = -2 x (side),  y = -2 x +  5 (side),
and make the assignments:

\(Clear[u, v, w, x, y, z]; 
    
{u[x_, y_, z_], v[x_, y_, z_], w[x_, y_, z_]} = {z - 2\ x, y - x, 
        y + 2\ x}\)

Note:
       z = 2  x -  1 ⟷ u[x, y, z] = -1,
       z = 2  x +  3 ⟷ u[x, y, z] = 3,
       y = x ⟷ v[x, y, z] = 0,
       y = x +  4 ⟷ v[x, y, z] = 4,
       y = -2 x ⟷ w[x, y, z] = 0,  and
       y = -2 x +  5 ⟷ w[x, y, z] = 5.
In uvw-space, (Subscript[R, uvw]) is everything inside and on the box bounded by the planes 
       u = -1, u = 3, v = 0, v = 4, w = 0, and w = 5.
Also
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\_xyz\)\%\ 
            \(e\^\(-x\)\) Sin[y] dx dy 
              dz\)\)\)
       \(\( = 
        ∫\_\%\(∫\_\%\(∫\_\(\ R\_uvw\)\%\ 
              \(e\^\(-x[u, v, w]\)\) Sin[y[u, v, w]] V\_xyz[u, v, w] 
                du dv 
                dw\)\)\)\)
        = 
        (∫_0^5)
            (∫_0^4)(∫_(-1)^3) 
              (e^(-x[u, v, w])) Sin[y[u, v, w]] (Subscript[V, xyz])[u, v, w] 
                du dv 
                dw.
To calculate x[u, v, w], y[u, v, w], z[u, v, w], and 
(Subscript[V, xyz])[u, v, w], use:

xyzsolved = 
      Solve[{u == u[x, y, z], v == v[x, y, z], w == w[x, y, z]}, {x, y, z}]

{x[u_, v_, w_], y[u_, v_, w_], z[u_, v_, w_]} = {x, y, z} /. 
        xyzsolved[[1]]

Now turn everything over to Mathematica:

\(Clear[gradx, grady, gradz, Vxyz]; 

    gradx[u_, v_, w_] = {D[x[u, v, w], u], D[x[u, v, w], v], 
        D[x[u, v, w], w]}; 

    grady[u_, v_, w_] = {D[y[u, v, w], u], D[y[u, v, w], v], 
        D[y[u, v, w], w]}; 

    gradz[u_, v_, w_] = {D[z[u, v, w], u], D[z[u, v, w], v], 
        D[z[u, v, w], w]}; 

    Vxyz[u_, v_, w_] = 
      Abs[Det[{gradx[u, v, w], grady[u, v, w], gradz[u, v, w]}]]\)

So
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\_xyz\)\%\ 
            \(e\^\(-x\)\) Sin[y] dx dy 
              dz\)\)\)
       \(\( = 
        ∫\_\%\(∫\_\%\(∫\_\(\ R\_uvw\)\%\ 
              \(e\^\(-x[u, v, w]\)\) Sin[y[u, v, w]] V\_xyz[u, v, w] 
                du dv 
                dw\)\)\)\)
        = 
        (∫_0^5)
            (∫_0^4)(∫_(-1)^3) 
              (e^(-x[u, v, w])) Sin[y[u, v, w]] ((1/3)) 
                du dv 
                dw.
The first integral 
       (∫_(-1)^3) 
        (e^(-x[u, v, w])) Sin[y[u, v, w]] ((1/3)) 
          du:

firstintegral = 
      (∫_(-1)^3) 
        (E^(-x[u, v, w])) Sin[y[u, v, w]] ((1/3)) \[DifferentialD]u

The second integral 
       (∫_0^4)(∫_(-1)^3) 
          (e^(-x[u, v, w])) Sin[y[u, v, w]] ((1/3)) du 
            dv is

secondintegral = (∫_0^4) firstintegral \[DifferentialD]v

The third and final integral
       (∫_0^5)
          (∫_0^4)(∫_(-1)^3) 
            (e^(-x[u, v, w])) Sin[y[u, v, w]] ((1/3)) 
              du dv dw 
is, in decimals,

NIntegrate[secondintegral, {w, 0, 5}]

You can get 
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\_xyz\)\%\ 
            \(e\^\(-x\)\) Sin[y] dx dy 
              dz\)\)\)
        = 
        (∫_0^5)
            (∫_0^4)(∫_(-1)^3) 
              (e^(-x[u, v, w])) Sin[y[u, v, w]] ((1/3)) 
                du dv 
                dw
in one step:

quickway = 
      NIntegrate[
        (1/3) (E^(-x[u, v, w])) Sin[y[u, v, w]], {w, 0, 5}, {v, 0, 4}, {u, 
          -1, 3}]

Went through that one like a hot knife going through butter.

B.2.b.iii)

The solid (Subscript[R, xyz]) consists of everything under the surface         
       z = 0.2  (e^(-y))
and over the region inside the ellipse 
       (x^2)+((((y/3)))^2)=1
in the xy-plane.
Here's a look:

\(Clear[x, y, z, r, t, s]; 
    
{x[r_, t_], y[r_, t_]} = {r\ Cos[t], 3\ r\ Sin[t]}; 

    \(z[r_, t_] = 0.2\ E\^\(-y[r, t]\); \)\){{rlow, rhigh}, {tlow, thigh}} = {{0, 1}, {0, 2 π}}; topskin = 
      ParametricPlot3D[{x[r, t], y[r, t], z[r, t]}, {r, rlow, rhigh}, {t, 
          tlow, thigh}, DisplayFunction -> Identity]; baseskin = 
      ParametricPlot3D[{x[r, t], y[r, t], 0}, {r, rlow, rhigh}, {t, tlow, 
          thigh}, PlotPoints -> {2, Automatic}, 
        DisplayFunction -> Identity]; \(sideskin = 
      ParametricPlot3D[{x[rhigh, t], y[rhigh, t], s\ z[rhigh, t]}, {s, 0, 
          1}, {t, tlow, thigh}, PlotPoints -> {2, Automatic}, 
        DisplayFunction -> Identity]; 


    Show[topskin, baseskin, sideskin, BoxRatios -> Automatic, 
      PlotRange -> All, AxesLabel -> {"x", "y", "z"}, 
      ViewPoint -> CMView, DisplayFunction -> $DisplayFunction]; 
    \)

Calculate
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\_xyz\)\%\ 
            \(e\^\(-y\)\) dx dy 
              dz\)\)\)
by transforming to an easy 3D integral.

Look at the plotting functions above:

{x[r, t], y[r, t]}

z[r, t]

Go to rst-space with:

\(Clear[x, y, z, r, s, t]; 
    
{x[r_, s_, t_], y[r_, s_, t_], z[r_, s_, t_]} = {r\ Cos[t], 
        3\ r\ Sin[t], 0.2\ s\ E\^\(-\((3\ r\ Sin[t])\)\)}\)

As you run 
       r from rlow to rhigh,
       t from tlow to thigh, and
       s from 0 to 1,
{x[r, s, t], y[r, s, t], z[r, s, t]} sweeps out all of (Subscript[R, xyz]); so (Subscript[R, rst])  is the box 
       rlow ≤ r ≤ rhigh, 0 ≤ s ≤ 1, tlow ≤ t ≤ thigh,
and
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\_xyz\)\%\ 
            \(e\^\(-y\)\) dx dy 
              dz\)\)\)
       \(\( = 
        ∫\_\%\(∫\_\%\(∫\_\(\ R\_rst\)\%\ 
              \(e\^\(-y[r, s, t]\)\) V\_xyz[r, s, t] dr 
                ds dt\)\)\)\)
        = 
        (∫_tlow^thigh)
            (∫_0^1)(∫_( rlow)^rhigh) 
              (e^(-y[r, s, t])) (Subscript[V, xyz])[r, s, t] dr 
                ds dt.

\(Clear[gradx, grady, gradz, Vxyz]; 

    gradx[r_, s_, t_] = {D[x[r, s, t], r], D[x[r, s, t], s], 
        D[x[r, s, t], t]}; 

    grady[r_, s_, t_] = {D[y[r, s, t], r], D[y[r, s, t], s], 
        D[y[r, s, t], t]}; 

    gradz[r_, s_, t_] = {D[z[r, s, t], r], D[z[r, s, t], s], 
        D[z[r, s, t], t]}; 

    Vxyz[r_, s_, t_] = 
      TrigReduce[Det[{gradx[r, s, t], grady[r, s, t], gradz[r, s, t]}]]\)

This is negative no matter what r, s, and t are in this box.  Reset it:

\(Clear[Vxyz]; 

    Vxyz[r_, s_, t_] = 
      \(-TrigReduce[Det[{gradx[r, s, t], grady[r, s, t], gradz[r, s, t]}]]
          \)\)

So
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\_xyz\)\%\ 
            \(e\^\(-y\)\) dx dy 
              dz\)\)\)
       \(\( = 
        ∫\_\%\(∫\_\%\(∫\_\(\ R\_rst\)\%\ 
              \(e\^\(-y[r, s, t]\)\) V\_xyz[r, s, t] dr 
                ds dt\)\)\)\)
        = 
        (∫_tlow^thigh)
            (∫_0^1)(∫_( rlow)^rhigh) 
              (e^(-y[r, s, t])) (Subscript[V, xyz])[r, s, t] dr 
                ds dt
is given by:

(∫_tlow^thigh)
        (∫_0^1)
          (∫_rlow^rhigh) (E^(-y[r, s, t])) Vxyz[r, s, t] 
            \[DifferentialD]r \[DifferentialD]s \[DifferentialD]t

Nice answer. 
But if you hadn't gone to rst-space to calculate the integral, getting this nice answer wouldn't have been nearly so simple.

B.3)  Mass and density

B.3.a)

Here x, y, and z are measured in centimeters.
The top of a solid R is the surface
       z = 32 -  (x^2) -  (y^2) 
and the bottom is the surface
       z = (x^2) +  (y^2) -  40.
The density of the solid varies from point to point, and is given by 
       density[x, y, z] = 0.3 ((x^2) +  (y^2) +  (z^2)) ((grams/(centimeter^3))) 
at {x, y, z}.  
It is fluffy near {0, 0, 0} and heavily dense far away from {0, 0, 0}.
Measure the total mass (in grams) of this solid.
If you dump it in a lake, will it sink or will it float?

The density function is what you integrate to calculate mass.  
As a result, the mass of a solid R whose density is density[x, y, z] (grams/(centimeter^3)) at a point {x, y, z} is measured by
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
            density[x, y, z] dx dy 
              dz\)\)\).

\(Clear[x, y, z, density]; 

    density[x_, y_, z_] = 0.03\ \((x\^2 + y\^2 + z\^2)\); \)

The top of R is the surface
       z = 32 -  (x^2) -  (y^2) 
and the bottom is the surface
       z = (x^2) +  (y^2) -  40.
With z so conveniently displayed, it is natural to integrate with respect to z first, so calculate
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
            density[x, y, z] dx dy 
              dz\)\)\)
       \(\( = 
        ∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
              density[x, y, z] dz dy 
                dx\)\)\)\)
To set the limits for the first integral,
       (∫_□^□) 
        density[x, y, z] dz,
fix x and y and enter the lowest and the highest values of z for this fixed x and y. The top skin of R is the surface
       z = 32 -  (x^2) -  (y^2) 
and the bottom is the surface
       z = (x^2) +  (y^2) -  40.

\(Clear[zhigh, zlow]; 
zhigh[x_, y_] = 32 - x\^2 - y\^2; 

    \(zlow[x_, y_] = x\^2 + y\^2 - 40; \)\)firstintegral = 
      (∫_(zlow[x, y])^(zhigh[x, y]))density[x, y, z] 
          \[DifferentialD]z

For the second integral, you have our choice whether to integrate with respect to y or to x.  
Go with y.
The main problem is deciding what limits to insert for yhigh[x] and ylow[x]  in the second integral
       (∫_(ylow[x])^(yhigh[x]
            ))(∫_(zlow[x, y])^(zhigh[x, y])) 
          density[x, y, z] dz dy.
The key to deciding how to set yhigh[x] and ylow[x] is to ask:
For a fixed x, what are the lowest and the highest values y can have no matter what z is?  
The simplest way of answering this is to plot the shadow of this solid on xy-paper with the y-axis running on the horizontal:

shadoweqn = Eliminate[{z == zhigh[x, y], z == zlow[x, y]}, z]

shadowcurves = Solve[shadoweqn, y]

\(Clear[ylow, yhigh]; 
    
{ylow[x_], yhigh[x_]} = {\(-\@\(36 - x\^2\)\), \@\(36 - x\^2\)}\)

shadowplot = 
      Plot[{ylow[x], yhigh[x]}, {x, -6, 6}, 
        PlotStyle -> {Blue, Red}, AxesLabel -> {"y", "x"}];

The second integral  
       (∫_(ylow[x])^(yhigh[x]
            ))(∫_(zlow[x, y])^(zhigh[x, y])) 
          density[x, y, z] dz dy
is:

secondintegral = 
      (∫_(ylow[x])^(yhigh[x]))Expand[firstintegral] 
          \[DifferentialD]y

For the third integral, you have no choice but to integrate with respect to x. 
The main problem is in deciding what limits to insert for xhigh and xlow in the third integral
       (∫_xlow^xhigh)
          (∫_(ylow[x])^(yhigh[x]
              ))(∫_(zlow[x, y])^(zhigh[x, y]))density[x, y, 
                z] dz dy 
              dx.

The key to deciding how to set xlow and xhigh is to ask:
What are the lowest and the highest values x can have no matter what z and y are?  
The simplest way of getting a hold on xlow and xhigh is to look at the shadow plot again:

Show[shadowplot];

In this one, xlow and xhigh appear where ylow[x] and yhigh[x] meet:

ends = Solve[yhigh[x] == ylow[x]]

That gives:

{xlow, xhigh} = {-6, 6};

The mass in grams is measured by the third integral
       (∫_xlow^xhigh)
          (∫_(ylow[x])^(yhigh[x]
              ))(∫_(zlow[x, y])^(zhigh[x, y])) 
            density[x, y, z] dz dy 
              dx:

mass = NIntegrate[secondintegral, {x, xlow, xhigh}]

To see whether this solid sinks or floats, calculate its volume:

volume = 
      (∫_xlow^xhigh)
          (∫_(ylow[x])^(yhigh[x]
              ))(∫_(zlow[x, y])^(zhigh[x, y]))1 
              \[DifferentialD]z \[DifferentialD]y \[DifferentialD]x

Compare:

{N[volume], mass}

One (centimeter^3) of water weighs one gram. 
So this solid outweighs the volume of water it displaces.
When you dump it in the lake, it sinks fast.

B.4)  Average value of a function

B.4.a.i) One variable

Given a function f[x] and an interval a ≤ x ≤ b, folks like to say that the average value of f[x] and this interval is given by 
       (1/((b  -  a)))(∫_a^b) f[x] dx.
Why do they like to say this?

They like to say this because it works.
To see it in action, go with 
       f[x] = Sin[x]  
and look at an experimental average of random values of f[x] for a = 0 ≤ x ≤ b = π:

\(Clear[f, x, k]; 
f[x_] = Sin[x]; 
\({a, b} = {0, π}; \)\)\(pointcount = 200; 

    experimentalaverage = 
      \(∑\+\(k = 1\)\%pointcount N[f[Random[Real, {N[a], N[b]}]]]
            \)\/pointcount\)

And the theoretical average value:

theoreticalaverage = 
      N[(((∫_a^b) f[x] \[DifferentialD]x)/(b - a))]

Compare:

{theoreticalaverage "= theoretical average", 
      experimentalaverage "= experimental average"}

See what happens when you use more points:

pointcount = 400; newexperimentalaverage = 
      (((∑_(k = 1)^pointcount) N[f[Random[Real, {N[a], N[b]}]]]
            )/pointcount); {theoreticalaverage "= theoretical average", 
      newexperimentalaverage "= experimental average"}

Uh-huh.
Reset the point count and rerun.

B.4.a.ii)

Wait a minute! This is supposed to be a mathematics course.  Where do you get off saying that something is correct because it works?

Good point.
You can explain why it works another way by looking at plots like this:

\(Clear[f, x]; 
f[x_] = Sin[x]; 
\({a, b} = {0, π}; \)\)fplot = 
      Plot[f[x], {x, a, b}, PlotStyle -> {{Thickness[0.01], Red}}, 
        DisplayFunction -> Identity]; average = N[(((∫_a^b) f[x] \[DifferentialD]x)/(b - a))]; 
    \(averagefplot = 
      Graphics[{Blue, 
          Polygon[{{a, 0}, {b, 0}, {b, average}, {a, average}}]}]; 


    Show[averagefplot, fplot, Axes -> True, 
      AxesLabel -> {"x", "y"}, PlotRange -> All, 
      AspectRatio -> Automatic, 
      DisplayFunction -> $DisplayFunction]; \)

The height of the rectangle is (1/((b  -  a)))(∫_a^b) f[x] dx.
The width of the rectangle is (b -  a).
The area inside the rectangle measures out to
       ((x] dx)/((b  -  a)))(b -  a)=(∫_a^b) f[x] dx
in square units.
This tells you that the area of the flat rectangle measures the (signed) area under the curve y = f[x]  and over the x-axis and for a ≤ x ≤ b.
This is why folks say that the average value of f[x] for a ≤ x ≤ b is the height of the rectangle, 
       (1/((b  -  a)))(∫_a^b) f[x] dx.
Try it again for
       f[x] = 2.3  (e^(-0.5 x)) for -1 = a ≤ x ≤ b = 2:

\(Clear[f, x]; 
f[x_] = 2.3\ E\^\(\(-0.5\)\ x\); 

    \({a, b} = {\(-1\), 2}; \)\)fplot = 
      Plot[f[x], {x, a, b}, PlotStyle -> {{Thickness[0.01], Red}}, 
        DisplayFunction -> Identity]; average = (((∫_a^b) f[x] \[DifferentialD]x)/(b - a)); 
    \(averagefplot = 
      Graphics[{Blue, 
          Polygon[{{a, 0}, {b, 0}, {b, average}, {a, average}}]}]; 


    Show[averagefplot, fplot, Axes -> True, 
      AxesLabel -> {"x", "y"}, PlotRange -> All, 
      AspectRatio -> Automatic, 
      DisplayFunction -> $DisplayFunction]; \)

And experimentally a couple of times:

pointcount = 50; experimentalaverage = 
      (((∑_(k = 1)^pointcount) N[f[Random[Real, {N[a], N[b]}]]]
            )/pointcount); {average "= theoretical average", 
      experimentalaverage "= experimental average"}

newpointcount = 300; newexperimentalaverage = 
      (((∑_(k = 1)^pointcount) N[f[Random[Real, {N[a], N[b]}]]]
            )/pointcount); {average "= theoretical average", 
      newexperimentalaverage "= experimental average"}

Usually not far off.

\(Clear[f, x]; 
f[x_] = E\^\(-x\); 
\({a, b} = {\(-1\), 2}; \)\)fplot = 
      Plot[f[x], {x, a, b}, PlotStyle -> {{Thickness[0.01], Red}}, 
        DisplayFunction -> Identity]; average = (((∫_a^b) f[x] \[DifferentialD]x)/(b - a)); 
    \(averagefplot = 
      Graphics[{Blue, 
          Polygon[{{a, 0}, {b, 0}, {b, average}, {a, average}}]}]; 


    Show[averagefplot, fplot, AxesLabel -> {"x", "y"}, 
      AspectRatio -> Automatic, 
      DisplayFunction -> $DisplayFunction]; \)

And experimentally several times:

theoreticalaverage = 
      N[(((∫_a^b) f[x] \[DifferentialD]x)/(b - a))]; pointcount = 200; experimentalaverage = 
      (((∑_(k = 1)^pointcount) N[f[Random[Real, {N[a], N[b]}]]]
            )/pointcount); {theoreticalaverage, experimentalaverage}

See what happens when you use more points:

pointcount = 200; experimentalaverage = 
      (((∑_(k = 1)^pointcount) N[f[Random[Real, {N[a], N[b]}]]]
            )/pointcount); {theoreticalaverage, experimentalaverage}

pointcount = 200; experimentalaverage = 
      (((∑_(k = 1)^pointcount) N[f[Random[Real, {N[a], N[b]}]]]
            )/pointcount); {theoreticalaverage, experimentalaverage}

pointcount = 200; experimentalaverage = 
      (((∑_(k = 1)^pointcount) N[f[Random[Real, {N[a], N[b]}]]]
            )/pointcount); {theoreticalaverage, experimentalaverage}

Usually not far off.

B.4.b) 2D averages.

Given a function f[x, y] and region R in xy-paper, folks like to say that the average value of f[x, y] on R is given by 
       (((∫_^)(∫_( R)^) 
                f[x, y] dx 
                  dy
                    )/((∫_^)(∫_( R
                  )^) dx dy))
Why do they like to say this?

Look at these:

{{a, b}, {c, d}} = {{0, 2}, {-1, 4}}; \(Clear[f, x, y]; 
f[x_, y_] = Sin[2\ x + y] + 1; 

    \(fplot = 
      Plot3D[f[x, y], {x, a, b}, {y, c, d}, Axes -> True, 
        AxesLabel -> {"x", "y", "z"}, 
        ViewPoint -> CMView, DisplayFunction -> Identity]; \)\)theoreticalaverage = 
      N[(((∫_a^b)
              (∫_c^d) f[x, y] \[DifferentialD]y 
                \[DifferentialD]x)/((∫_a^b)
              (∫_c^d) 1 \[DifferentialD]y \[DifferentialD]x))]; 
    \(averagefplot = 
      Plot3D[theoreticalaverage, {x, a, b}, {y, c, d}, 
        PlotPoints -> {2, 2}, Axes -> True, 
        AxesLabel -> {"x", "y", "z"}, 
        ViewPoint -> CMView, DisplayFunction -> Identity]; 


    Show[GraphicsArray[{fplot, averagefplot}], 
      DisplayFunction -> $DisplayFunction]; \)

The volume measurements of the solids between the xy-plane and the plotted surfaces are the same.
Reason:
       volume left = 
        (∫_^)(∫_( R)^) 
            f[x, y] dx dy
       volume right = base times height
       =(∫_^)(∫_( R)^) dxdy(\(\(∫\_\%\(∫\_R\%f
                            \([x, y]\) dx 
                            dy
                              \)\)\/\(∫\_\%\(∫\_R\%dx dy\)\)\))
        = 
        (∫_^)(∫_( R)^) f[x, y] 
              dx dy.
Both solids have the same base, so the height of the solid on the right must be the average value of f[x, y] for {x, y} in R.
See what happens experimentally:

pointcount = 300; experimentalaverage = 
      (1/pointcount) 
        ((∑_(k = 1)^pointcount) N[
              f[Random[Real, {N[a], N[b]}], Random[Real, {N[c], N[d]}]]]); 
    {theoreticalaverage "= theoretical average", 
      experimentalaverage "= experimental average"}

Again:

pointcount = 500; experimentalaverage = 
      (1/pointcount) 
        ((∑_(k = 1)^pointcount) N[
              f[Random[Real, {N[a], N[b]}], Random[Real, {N[c], N[d]}]]]); 
    {theoreticalaverage "= theoretical average", 
      experimentalaverage "= experimental average"}

Lookin' good and feeling good.

B.4.c) 3D averages.

Given a function f[x, y, z] and a solid R in xyz-space, folks like to say that the average value of f[x, y, z] on R is given by 
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
            f[x, y, z] dx dy 
              dz/
                \(∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
                      dx dy 
                        dz\)\)\)\)\)\)
Why do they like to say this?

This time, a plot is not possible.
Try an experiment with R as the box
       0 ≤ x ≤ 2, -1 ≤ y ≤ 4, and -2 ≤ z ≤ 3:

{{xlow, xhigh}, {ylow, yhigh}, {zlow, zhigh}} = {{0, 2}, {-1, 
          4}, {-2, 3}}; \(Clear[f, x, y, z]; 
f[x_, y_, z_] = Sin[2\ x + y - z] + 1; 

    \(theoreticalaverage = 
      N[\(∫\_xlow\%xhigh
              \(∫\_ylow\%yhigh
                \(∫\_zlow\%zhigh f[x, y, z] \[DifferentialD]z 
                  \[DifferentialD]y 
                  \[DifferentialD]x\)\)\)\/\(∫\_xlow\%xhigh
              \(∫\_ylow\%yhigh
                \(∫\_zlow\%zhigh 1 \[DifferentialD]z 
                  \[DifferentialD]y \[DifferentialD]x\)\)\)]; \)\)pointcount = 100; experimentalaverage = 
      (1/pointcount) 
        ((∑_(k = 1)^pointcount) N[
              f[Random[Real, {N[xlow], N[xhigh]}], 
                Random[Real, {N[ylow], N[yhigh]}], 
                Random[Real, {N[zlow], N[zhigh]}]]]); {theoreticalaverage "= theoretical average", 
      experimentalaverage "= experimentalaverage"}

pointcount = 200; newexperimentalaverage = 
      (1/pointcount) 
        ((∑_(k = 1)^pointcount) N[
              f[Random[Real, {N[xlow], N[xhigh]}], 
                Random[Real, {N[ylow], N[yhigh]}], 
                Random[Real, {N[zlow], N[zhigh]}]]]); {theoreticalaverage "= theoretical average", 
      newexperimentalaverage "= experimentalaverage"}

And look at this:

{theoreticalaverage "= theoretical average", 
      (1/2) (experimentalaverage + newexperimentalaverage) 
        "= averageexperimentalaverage"}

Nice.
Another way to look at this is to note that when you go with
       \(average = 
        ∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
              f[x, y, z] dx dy 
                dz/
                  \(∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
                        dx dy 
                          dz\)\)\)\)\)\), 
then you get
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
              f[x, y, z] dx dy 
                dz\)\) = 
        ∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
              average dx dy 
                dz\)\)\).