2.05 2D Integrals and the Gauss-Green Formula
 BASICS

Off[General::spell]; Off[General::spell1]; Off[Plot::plnr]; Off[ParametricPlot3D::ppcom]; SetOptions[Limit, Analytic -> True]; \(Needs["Graphics`FilledPlot`"]; 


    If[MemberQ[{"Global`", \ "Miscellaneous`SIUnits`", \ 
          "Graphics`Colors`"}, 
\ \ \ \ Context[Gray]], \ 
      Remove["*`Gray"]]; 
<< "Miscellaneous`Units`"; 

    Remove["*`Gray"]; 
<< "Graphics`Colors`"; 

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

    
$Pre\  = \ Chop; 
CMView\  = \ {2.7, \ 1.6, \ 1.2}; \)

B.1)  2D integrals ∫∫(Subscript[, R]) f[x, y]dxdy  for volume measurements

Here's a plot of a function y = f[x]:

\(Clear[f, x]; 
f[x_] = E\^\(\(-0.2\)\ x\)\ Cos[4\ x]; 
    
{a, b} = {0, 4}; 


    fplot = Plot[f[x], {x, a, b}, 
        PlotStyle -> {{Blue, Thickness[0.01]}}, 
        AxesLabel -> {"x", "y"}]; \)

To get the points on the plot, you go to a point 
       {(Subscript[x, 0]), 0} 
on the x-axis underneath or above the plot and then you run a line to the point
       {(Subscript[x, 0]), f[(Subscript[x, 0])]} 
on the curve like this:

\(xo = 1.5; 


    Show[fplot, 
	Graphics[{PointSize[0.02], Point[{xo, 0}]}], 
      Graphics[{PointSize[0.02], Point[{xo, f[xo]}]}], 
      Graphics[{Red, Line[{{xo, 0}, {xo, f[xo]}}]}]]; \)

Play with this by resetting (Subscript[x, 0]) and rerunning.
You can get a pretty good idea of what the curve looks like by just looking at the tips of the little line segments:

\(xjump = \(b - a\)\/50; 


    Show[fplot, 
	
      Table[{Graphics[{PointSize[0.02], Point[{x, 0}]}], 
          Graphics[{PointSize[0.02], Point[{x, f[x]}]}], 
          Graphics[{Red, Line[{{x, 0}, {x, f[x]}}]}]}, {x, a, b, xjump}]]; 
    \)

You can plot a surface z = f[x, y] like this:

\(\(Clear[f, x, y]; 
f[x_, y_]\  = \ 3.1\ x\^2 + 2.3\ y\^2; 
    
{{a, b}, {c, d}}\  = \ {{\(-2\), 3}, {\(-1\), 4}}; 


    surfaceplot\  = \ 
      Plot3D[f[x, y], {x, a, b}, {y, c, d}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        DisplayFunction -> Identity]; 
\ \ \ 
spacer\  = \ 0.2; 

    threedims\  = \ 
\ \ \ \ \ \ 
      Graphics3D[{

          \ \ \ \ \ \ \ \ \ \ \ \ {Blue, Line[{{a, 0, 0}, {b, 0, 0}}]}, 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          Text["x", {b\  + \ spacer, 0, 0}], 

          \ \ \ \ \ \ \ \ \ \ \ \ {Blue, Line[{{0, c, 0}, {0, d, 0}}]}, 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          Text["y", {0, \ d\  + \ spacer, 0}], 

          \ \ \ \ \ \ \ \ \ \ \ \ {Blue, Line[{{0, 0, 0}, {0, 0, 50}}]}, 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          Text["z", {0, 0, 50\  + \ spacer}]}]; 

    \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
CMView\  = \ {2.7, 1.6, 1.2}; 

    fplot\  = \ 
      Show[surfaceplot, threedims, 
\ \ \ \ \ \ \ \ ViewPoint -> CMView, 
        PlotRange -> All, 
\ \ \ \ \ \ \ \ 
        DisplayFunction -> $DisplayFunction]; \ 
\)\)

B.1.a)

Explain the meaning of the plotted points that make up the surface.

You do it just as you did it above.
To get the points on the plot, you go to a point
       {(Subscript[x, 0]), (Subscript[y, 0]), 0} 
on the xy-plane underneath or above the surface. Then you run a line to the 3 D point
       {(Subscript[x, 0]), (Subscript[y, 0]), f[(Subscript[x, 0]), (Subscript[y, 0])]} 
like this:

\({xo, yo} = {2.5, 3}; 


    Show[fplot, 
	Graphics3D[{PointSize[0.02], Point[{xo, yo, 0}]}], 
	
      Graphics3D[{PointSize[0.02], Point[{xo, yo, f[xo, yo]}]}], 
	
      Graphics3D[{Red, Line[{{xo, yo, 0}, {xo, yo, f[xo, yo]}}]}]]; \)

Play with this by resetting (Subscript[x, 0]) and (Subscript[y, 0]) and rerunning.
You can get a pretty good idea of what the surface looks like by just looking at the tips of the little line segments:

xjump = ((b - a)/8); \(yjump = \(d - c\)\/8; 


    Show[fplot, 
	
      Table[{Graphics3D[{PointSize[0.02], Point[{x, y, 0}]}], 
			
          Graphics3D[{PointSize[0.02], Point[{x, y, f[x, y]}]}], 
			
          Graphics3D[{Red, Line[{{x, y, 0}, {x, y, f[x, y]}}]}]}, {x, a, b, 
          xjump}, {y, c, d, yjump}]]; \)

Think of the line segments as poles propping up the surface.
If the base of the pole is at {x, y}, then the height of the pole is f[x, y].

B.1.b.i) How to use a double integral to measure volume

Take another look at the plot from part a) immediately above:

\(\(Clear[f, x, y]; 
f[x_, y_]\  = \ 3.1\ x\^2 + 2.3\ y\^2; 
    
{{a, b}, {c, d}}\  = \ {{\(-2\), 3}, {\(-1\), 4}}; 


    surfaceplot\  = \ 
\ \ \ 
      Plot3D[f[x, y], {x, a, b}, {y, c, d}, 
\ \ \ 
        DisplayFunction -> Identity]; 
\ \ \ 
spacer\  = \ 0.2; 

    threedims\  = \ 
\ \ 
      Graphics3D[{
\ \ \ \ \ {Blue, Line[{{a, 0, 0}, {b, 0, 0}}]}, 

          \ \ \ \ \ \ \ \ \ \ \ \ \ Text["x", {b\  + \ spacer, 0, 0}], 

          \ \ \ \ \ {Blue, Line[{{0, c, 0}, {0, d, 0}}]}, 

          \ \ \ \ \ \ \ \ \ \ \ \ \ Text["y", {0, \ d\  + \ spacer, 0}], 
          
\ \ \ \ \ {Blue, Line[{{0, 0, 0}, {0, 0, 50}}]}, 

          \ \ \ \ \ \ \ \ \ \ \ \ \ Text["z", {0, 0, 50\  + \ spacer}]}]; 
    
CMView\  = \ {2.7, 1.6, 1.2}; 


    fplot\  = \ 
\ \ \ 
      Show[surfaceplot, threedims, 
\ \ \ ViewPoint -> CMView, 
        PlotRange -> All, 
\ \ \ DisplayFunction -> $DisplayFunction]; \ 
    \)\)

How do you use the double integral
       ∫(∫_R) f[x, y] dx 
            dy 
to measure the volume of the solid consisting of everything between the surface plotted above and the xy-plane?

Look again at the poles like those plotted above:

xjump = ((b - a)/10); \(yjump = \(d - c\)\/10; 


    Show[fplot, 
	
      Table[{Graphics3D[{PointSize[0.02], Point[{x, y, 0}]}], 
			
          Graphics3D[{Red, Line[{{x, y, 0}, {x, y, f[x, y]}}]}]}, {x, a, b, 
          xjump}, {y, c, d, yjump}]]; \)

The solid region you're talking about consists of all the poles you get when you plot all the possible poles (not just those plotted above).
When you plotted this surface on the top, you ran 
       x from a to b
and you ran
       y from c to d.  
The volume of the solid is given by the double integral
       ∫(∫_R) f[x, y] dx 
            dy 
where R is the rectangle in the xy-plane consisting of all points {x, y, 0} with a ≤ x ≤ b and c ≤ y ≤ d.
You can calculate the volume measurement
       ∫(∫_R) f[x, y] dx 
            dy 
by writing it as
       (∫_c^d)
          (∫_a^b) f[x, y] dx 
            dy
which you get from Mathematica with:

(∫_a^b)
        (∫_c^d) f[x, y] \[DifferentialD]y \[DifferentialD]x

The volume of the solid measures out to 430 cubic units.

B.1.b.ii) Why double integrals measure volume

What is the physical process that confirms that the double integral
       ∫(∫_R) f[x, y] dx 
              dy = 
        (∫_c^d)
            (∫_a^b) f[x, y] dx 
              dy
actually measures the volume of the solid described in part i)?

Look at some slices of the solid consisting of everything between the surface plotted above and the xy-plane:

\(Clear[slice, x, y, t]; 

    slice[y_] := 
      ParametricPlot3D[{x, y, 0} + t\ {0, 0, f[x, y]}, {x, a, b}, {t, 0, 1}, 
        PlotPoints -> {Automatic, 2}, DisplayFunction -> Identity]; 
    
jump = \(d - c\)\/3; 


    Table[Show[fplot, slice[y], DisplayFunction -> $DisplayFunction], {
        y, c, d, jump}]; \)

The physical process of calculating 
       (∫_c^d)
          (∫_a^b) f[x, y] dx 
            dy 
involves two steps.
First step:
For each y with c ≤ y ≤ d, you calculate
       (∫_a^b) f[x, y] dx:

\(Clear[first, x, y]; 

    first[y_] = ∫\_a\%b f[x, y] \[DifferentialD]x\)

As y advances from c to d, first[y] measures the area of the slices that sweep out the solid. 
For instance,

first[3]

measures the area of this slice of the solid:

Show[fplot, slice[3], DisplayFunction -> $DisplayFunction];

Second step: 
Calculate (∫_c^d) first[y] dy = 
        (∫_c^d)
            (∫_a^b) f[x, y] dx 
              dy.
This second integral sums up the area measurements of the individual slices and measures the volume.

second = (∫_c^d) first[y] \[DifferentialD]y

Compare to Mathematica's calculation of
       (∫_c^d)
          (∫_a^b) f[x, y] dx 
            dy:

(∫_c^d)
        (∫_a^b) f[x, y] \[DifferentialD]x \[DifferentialD]y

Bingo.

B.2)  What ∫(∫_R) f[x, y] dx 
            dy means when f[x, y] isn't always positive

B.2.a)

Here's a calculation of 
       (∫_0^4)
          (∫_0^2) f[x, y] dx 
            dy 
for f[x, y] = (x - 3) (E^(0.43  y)):

\(Clear[f, x, y]; 
f[x_, y_] = \((x - 3)\)\ E\^\(0.43\ y\); 

    ∫\_0\%2
        \(∫\_0\%4 f[x, y] \[DifferentialD]y \[DifferentialD]x\)\)

Negative.
When the dweeby Calculus Cal looked at this, he said, "This has got to be wrong because double integrals measure volume and volume can't be negative."
What do you say to Cal in response?

First, you tell Cal to stop picking his nose.
But then you realize that Cal's question does have scientific merit and you answer it by telling Cal to look at the integral
       (∫_0^4)
          (∫_0^2) f[x, y] dx 
            dy  
and the corresponding plot of f[x, y]:

\({a, b}\  = {0, 2}; 
{c, d}\  = \ {0, 4}; 

    fplot\  = \ 
      Plot3D[f[x, y], {x, a, b}, {y, c, d}, 
\ \ \ \ \ \ \ \ 
        DisplayFunction -> Identity]; 

spacer\  = \ 0.2; 

    threedims\  = \ 
\ \ \ \ \ \ 
      Graphics3D[{

          \ \ \ \ \ \ \ \ \ \ \ \ {Blue, Line[{{a, 0, 0}, {b, 0, 0}}]}, 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          Text["x", {b\  + \ spacer, 0, 0}], 

          \ \ \ \ \ \ \ \ \ \ \ \ {Blue, Line[{{0, c, 0}, {0, d, 0}}]}, 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          Text["y", {0, \ d\  + \ spacer, 0}], 

          \ \ \ \ \ \ \ \ \ \ \ \ {Blue, Line[{{0, 0, \(-10\)}, {0, 0, 5}}]}, 
          
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          Text["z", {0, 0, 5\  + \ spacer}]}]; 

    CMView\  = \ {2.7, 1.6, 1.2}; 


    Show[fplot, threedims, ViewPoint -> CMView, 
\ \ \ \ \ 
      DisplayFunction -> $DisplayFunction]; \)

You say that the plot reveals that f[x, y] is negative for the x's and y's involved in the calculation of
       (∫_0^4)
          (∫_0^2) f[x, y] dx 
            dy. 
This is why (∫_0^4)
          (∫_0^2) f[x, y] dx 
            dy is negative.
You go on to tell Cal that if he wants a measurement of the volume of the solid whose BOTTOM skin is the plotted part of the surface and whose TOP skin is the rectangle in the xy-plane consisting of the points {x, y, 0} with 0 ≤ x ≤ 2 and 0 ≤ y ≤ 4, then he can get the measurement by calculating
       -
        (∫_0^4)
            (∫_0^2) f[x, y] dx 
              dy:

-(∫_0^2)
          (∫_0^4) f[x, y] \[DifferentialD]y 
            \[DifferentialD]x

Positive.
Now Cal is happy and so are you.

B.2.b.i)

Here's a calculation of
       (∫_0^4)
          (∫_0^5) f[x, y] dx 
            dy 
for f[x, y] = (x - 3) (E^(0.43  y)):

\(Clear[f, x, y]; 
f[x_, y_]\  = \((x - 3)\)\ E\^\(0.43\ y\); 

    Integrate[f[x, y], {x, 0, 5}, {y, 0, 4}]\)

Try to interpret this calculation.

Look at the integral
       (∫_0^4)
          (∫_0^5) f[x, y] dx 
            dy   
and the corresponding plot of f[x, y]:

\({a, b}\  = {0, 5}; 
{c, d}\  = \ {0, 4}; 

    fplot\  = \ 
      Plot3D[f[x, y], {x, a, b}, {y, c, d}, 
\ \ \ \ \ \ \ \ 
        DisplayFunction -> Identity]; 

spacer\  = \ 0.2; 

    threedims\  = \ 
\ \ \ \ \ \ 
      Graphics3D[{

          \ \ \ \ \ \ \ \ \ \ \ \ {Blue, Line[{{a, 0, 0}, {b, 0, 0}}]}, 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          Text["x", {b\  + \ spacer, 0, 0}], 

          \ \ \ \ \ \ \ \ \ \ \ \ {Blue, Line[{{0, c, 0}, {0, d, 0}}]}, 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          Text["y", {0, \ d\  + \ spacer, 0}], 

          \ \ \ \ \ \ \ \ \ \ \ \ {Blue, Line[{{0, 0, \(-10\)}, {0, 0, 5}}]}, 
          
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          Text["z", {0, 0, 5\  + \ spacer}]}]; 

    CMView\  = \ {2.7, 1.6, 1.2}; 


    Show[fplot, threedims, ViewPoint -> CMView, 
\ \ \ \ \ 
      DisplayFunction -> $DisplayFunction]; \)

Look at the plot to see that f[x, y] is sometimes positive and sometimes negative as x varies between 0 and 5 and y varies between 0 and 4.
The upshot:
       (∫_0^4)
          (∫_0^5) f[x, y] dx 
            dy   
does not measure any specific volume.

B.2.b.ii)

Now you know that when you go with
       f[x, y] = (x - 3) (E^(0.43  y)),
 then you can calculate
       (∫_0^4)
          (∫_0^5) f[x, y] dx 
            dy.

\(Clear[f, x, y]; 
f[x_, y_] = \((x - 3)\)\ E\^\(0.43\ y\); 


    ∫\_0\%5
        \(∫\_0\%4 f[x, y] \[DifferentialD]y \[DifferentialD]x\)\)

But the integral doesn't measure any specific volume because f[x, y] takes on both positive and negative values as x varies between 0 and 5 and y varies between 0 and 4.
Just how can you interpret this integral?

Look at:

f[x, y]

You can see that
       f[x, y] ≥ 0 for x ≥ 3 
and
       f[x, y] ≤ 0 for x ≤ 3.
To get the interpretation of (∫_0^4)
          (∫_0^5) f[x, y] dx 
            dy, write
       (∫_0^4)
            (∫_0^5) f[x, y] dx 
              dy = 
       (∫_0^4)
            (∫_0^3) f[x, y] dx 
              dy + 
        (∫_0^4)
            (∫_3^5) f[x, y] dx 
              dy.
→   In the first integral, f[x, y] ≥ 0; so
       (∫_0^4)
          (∫_0^3) f[x, y] dx 
            dy 
measures the volume of the solid whose BOTTOM skin is the rectangle in the xy-plane consisting of all points {x, y, 0} with 0 ≤ x ≤ 3 and 0 ≤ y ≤ 4 and whose TOP skin is the surface z = f[x, y] plotted with 0 ≤ x ≤ 3 and 0 ≤ y ≤ 4.
→   In the second integral, f[x, y] ≤ 0; so
       (∫_0^4)
          (∫_3^5) f[x, y] dx 
            dy 
measures the NEGATIVE of the volume of the solid whose TOP skin is the rectangle in the xy-plane consisting of all points {x, y, 0} with 3 ≤ x ≤ 5 and 0 ≤ y ≤ 4 and whose BOTTOM skin is the surface z = f[x, y] plotted with 3 ≤ x ≤ 5 and 0 ≤ y ≤ 4.
Consequently,
       (∫_0^4)
            (∫_0^5) f[x, y] dx 
              dy = 
        (∫_0^4)
              (∫_0^3) f[x, y] dx 
                dy + 
          (∫_0^4)
              (∫_3^5) f[x, y] dx 
                dy
measures the difference of the two volumes of the two solids calculated above.
Just a little bit esoteric.

B.3)  Trying to calculate ∫(∫_R) f[x, y] dx 
            dy when R isn't a rectangle

B.3.a) Low to high

Here is a little region R on the xy-plane.

\(Clear[high, low, x]; 
high[x_] = 2\ Sin[\(π\ x\)\/3]; 

    low[x_] = 1\/2\ x\ \((x - 3)\); 


    regionplot = 
      Plot[{high[x], low[x]}, {x, 0, 3}, 
        PlotStyle -> {{Blue, Thickness[0.01]}}, 
        AspectRatio -> Automatic, AxesLabel -> {"x", "y"}, 
        Epilog -> {Text["y = high[x]", {0.7, 1.5}], 
            Text["y = low[x]", {2.0, \(-0.8\)}], 
            Text["R", {1.5, 1}]}]; \)

The region R consists of all points between and on the plotted curves. Calculate
       ∫(∫_R) f[x, y] dx 
            dy
for
       f[x, y] = (E^(-0.4 x  -  0.3  y)).

Enter the function f[x, y]:

\(Clear[f, x, y]; 
f[x_, y_] = E\^\(\(-0.4\)\ x - 0.3\ y\)\)

To set up 
       ∫(∫_R) f[x, y] dx 
            dy
for calculation, you need a convenient way of slicing R.
Here's one:

\(Clear[slice, sliceplot, x]; 

    slice[x_] := 
      Graphics[{Red, Thickness[0.01], Line[{{x, low[x]}, {x, high[x]}}]}]; 

    sliceplot[x_] := 
      Show[regionplot, slice[x], DisplayFunction -> Identity]; 


    Show[sliceplot[1], DisplayFunction -> $DisplayFunction]; \)

You slice from low to high, and this means that you integrate with respect to y first. This tells you how to set up
       ∫(∫_R) f[x, y] dx 
              dy = 
        (∫_0^3)
            ((∫_(low[x])^(high[x]))f[x, y] dy)
               dx 
for calculation:

\(Clear[first]; 

    first[x_] = 
      ∫\_\(low[x]\)\%\(high[x]\)f[x, y] \[DifferentialD]y\)

second = NIntegrate[first[x], {x, 0, 3}]

Here's how to get Mathematica to calculate the double integral
       ∫(∫_R) f[x, y] dx 
              dy = 
        (∫_0^3)
            (∫_(low[x])^(high[x]))f[x, y] dy 
              dx
in one sweet instruction with NIntegrate:

NIntegrate[f[x, y], {x, 0, 3}, {y, low[x], high[x]}]

Done.

B.3.b) Left to right

Here is another little region on the xy-plane.

\(Clear[left, right, y]; 
left[y_] = y; 
right[y_] = 2 - y\^2; 


    regionplot = 
      ParametricPlot[{{left[y], y}, {right[y], y}}, {y, \(-2\), 1}, 
        PlotStyle -> {{Blue, Thickness[0.01]}}, 
        AspectRatio -> Automatic, AxesLabel -> {"x", "y"}, 
        Epilog -> {Text["x = left[y]", {\(-1\), \(-0.6\)}], 
            Text["x = right[y]", {1.5, \(-1\)}], 
            Text["R", {1, \(-0.5\)}]}]; \)

The region R consists of all points between and on the plotted curves. Calculate
       ∫(∫_R) f[x, y] dx 
            dy 
for
       f[x, y] = x y.

Enter the function f[x, y]:

\(Clear[f, x, y]; 
f[x_, y_] = x\ y\)

To set up
       ∫(∫_R) f[x, y] dx 
            dy 
for calculation, you need a convenient way of slicing R.
Here's one:

\(Clear[slice, sliceplot, y]; 

    slice[y_] := 
      Graphics[{Red, Thickness[0.01], Line[{{left[y], y}, {right[y], y}}]}]; 
    
sliceplot[y_] := 
      Show[regionplot, slice[y], DisplayFunction -> Identity]; 


    Show[sliceplot[\(-0.5\)], DisplayFunction -> $DisplayFunction]; \)

You slice from left to right, and this means that you integrate with respect to x first. This tells you how to set up
       ∫(∫_R) f[x, y] dx 
              dy = 
        (∫_(-2)^1)
            (∫_(left[y])^(right[y]))f[x, y] dx 
              dy 
for calculation:

\(Clear[first]; 

    first[y_] = 
      ∫\_\(left[y]\)\%\(right[y]\)f[x, y] \[DifferentialD]x\)

second = (∫_(-2)^1) first[x] \[DifferentialD]x

Here's how to get Mathematica to calculate the double integral
       ∫(∫_R) f[x, y] dx 
              dy = 
        (∫_(-2)^1)
            (∫_(left[y])^(right[y]))f[x, y] dx 
              dy 
in one sweet instruction:

(∫_(-2)^1)
        (∫_(left[y])^(right[y]))f[x, y] \[DifferentialD]x 
          \[DifferentialD]y

Done.

B.4)  The Gauss-Green formula helps you calculate
        ∫(∫_R) f[x, y] dx 
            dy

One part of the art of mathematics is knowing how to avoid miserable calculations.  The Gauss-Green formula is a tool for avoiding miserable calculations.
The Gauss-Green formula says this: 
When you have a parameterization {x[t], y[t]} of the boundary of a region R in the xy-plane with the extra feature that as t advances from a to b, the points {x[t], y[t]} advance around the boundary of R exactly once in the counterclockwise way, then you are guaranteed that
       ∫(∫_R)(∂_x) n[x, y] - 
        (∂_y) m[x, y] dx dy
        = 
        (∫_a^b) m[x[t], y[t]] x'[t] +  
          n[x[t], y[t]] y'[t] dt.
It's a funny-looking formula, but don't let it intimidate you. Its weirdness gives it the utility of a crescent wrench because you can adjust it to many situations. 
Get to feel comfortable with this formula because you can use it to side-step a lot of miserable calculations.

B.4.a) Gauss-Green for calculating double integrals

Use the Gauss-Green formula to help calculate the double integral
       ∫(∫_R)((x^2) + (y^2)) (E^(-y)) 
            dx dy
where R is the region inside the ellipse
       ((((x/2)))^2)+(y^2)=1.

Here's how you use the Gauss-Green formula
       ∫(∫_R)(∂_x) n[x, y] - 
        (∂_y) m[x, y] dx dy
        = 
        (∫_a^b) m[x[t], y[t]] x'[t] + 
          n[x[t], y[t]] y'[t] dt 
to calculate
       ∫(∫_R)((x^2) + (y^2)) (E^(-y)) 
            dx dy.
You just say
       f[x, y] = ((x^2) + (y^2)) (E^(-y))
       m[x, y] = 0 and
       n[x, y] = (∫_0^x) f[s, y] ds:

\(Clear[x, y, n, m, t, s]; 
f[x_, y_] = \((x\^2 + y\^2)\)\ E\^\(-y\); 

    \(m[x_, y_] = 0; \)\)n[x_, y_] = (∫_0^x) f[s, y] \[DifferentialD]s

This gives you
       f[x, y] = D[n[x, y], x] - D[m[x, y], y]:

{D[n[x, y], x] -  D[m[x, y], y],  f[x, y]}

Now you know that
       ∫(∫_R)((x^2) + (y^2)) (E^(-y)) 
              dx dy = 
        ∫(∫_R) f[x, y] dx 
              dy
        = 
        ∫(∫_R)(D[n[x, y], x] - D[m[x, y], y]) 
              dx dy
        = 
        (∫_a^b) m[x[t], y[t]] x'[t] + 
          n[x[t], y[t]] y'[t] dt 
where {x[t], y[t]} is a parameterization of the ellipse
       ((((x/2)))^2)+(y^2)=1 
chosen so that the points {x[t], y[t]} advance around the ellipse exactly once in the counterclockwise fashion as t advances from a to b.
One parameterization of
       ((((x/2)))^2)+(y^2)=1 
that does this nicely is
       {x[t], y[t]} = {2  Cos[t], Sin[t]}
with a = 0 and b = 2  π.
The upshot:
       ∫(∫_R)((x^2) + (y^2)) (E^(-y)) 
            dx dy
is given by:

\(a\  = \ 0; 
b\  = \ 2\ Pi; 

Clear[x, y, t]; 
    
{x[t_], y[t_]}\  = \ {2\ Cos[t], Sin[t]}; 


    NIntegrate[

      m[x[t], y[t]]\ \(x'\)[t]\  + \ n[x[t], y[t]]\ \(y'\)[t], {t, a, b}]\)

That's all there is to it. 
You're out of here without worrying about slicing or anything like that.

B.4.b)

To use the Gauss-Green formula to calculate
       ∫(∫_R) f[x, y] dx 
            dy,  
you need a parameterization {x[t], y[t]} of the boundary R with the extra feature that as t advances from a to b, the points {x[t], y[t]} advance around the boundary of R exactly once in the counterclockwise way. Then you are guaranteed that when you set
       n[x, y] = (∫_0^x) f[s, y] ds, 
you are guaranteed that
       ∫(∫_R) f[x, y] dx 
              dy = 
        (∫_a^b) n[x[t], y[t]] y'[t] dt.
The beauty of this is that it allows you to calculate a 2 D integral by means of a single integral that you can usually calculate with NIntegrate.
Discuss why this works.

It's all a matter of the fundamental formula of calculus.  
The fundamental formula of calculus tells you that when you start with a function f[x] and you put
       n[x] = (∫_0^x) f[s] ds, 
you are guaranteed that
       D[n[x], x] = f[x].
Try it out.

\(Clear[f, n, s, x]; 
f[x_] = x\ Sin[Random[Real, {\(-5\), 5}]\ x]; 


    n[x_] = ∫\_0\%x f[s] \[DifferentialD]s\)

Compare D[n[x], x] and f[x]:

{D[n[x], x], f[x]}

The same. Play with this by changing f[x] and rerunning.
The fundamental formula also says that when you start with a function f[x, y] and you put
       n[x, y] = (∫_0^x) f[s, y] ds
you are guaranteed that
       (∂_y) n[x, y] = f[x, y].
Try it out.

\(Clear[f, n, s, x, y]; 
f[x_, y_] = x\^2\ Sin[y]; 


    n[x_, y_] = ∫\_0\%x f[s, y] \[DifferentialD]s\)

Compare D[n[x, y], x] and f[x, y]:

{D[n[x, y], x], f[x, y]}

The same. Play with this by changing f[x, y] and rerunning.
When you feel comfortable with this, then you are ready to use the Gauss-Green formula comfortably.

Here is the Gauss-Green formula:
When you have a parameterization {x[t], y[t]} of the boundary of a region R in the xy-plane with the extra feature that as t advances from a to b, the points {x[t], y[t]} advance around the boundary of R exactly once in the counterclockwise way, then you are guaranteed that
       ∫(∫_R)(∂_x) n[x, y] - 
        (∂_y) m[x, y] dx dy
        = 
        (∫_a^b) m[x[t], y[t]] x'[t] + 
          n[x[t], y[t]] y'[t] dt.
To use this formula to calculate
       ∫(∫_R) f[x, y] dx 
            dy 
for a given function f[x, y], you just set
       m[x, y] = 0 for all x's and y's
and you set
       n[x, y] = (∫_0^x) f[s, y] ds.
This guarantees that
      (∂_x) n[x, y] = f[x, y] and (∂_y) m[x, y] = 0.
So
       ∫(∫_R) f[x, y] dx 
            dy
        = 
        ∫(∫_R) D[n[x, y], x] dx 
              dy
        = 
        ∫(∫_R)(D[n[x, y], x] - 0) dx 
              dy
        = 
        ∫(∫_R)(D[n[x, y], x] - D[m[x, y], y]) 
              dx dy
        = 
        (∫_a^b) m[x[t], y[t]] x'[t] + 
          n[x[t], y[t]] y'[t] dt
        = 
        (∫_a^b)(0 + n[x[t], y[t]] y'[t]) 
            dt
        = 
        (∫_a^b) n[x[t], y[t]] y'[t] dt.
Explanation complete.

B.4.c)

The region R in the xy-plane consists of everything inside and on the curve you see below:

\(Clear[x, y, t]; 
    
{x[t_], 
        y[t_]} = {2, 3} + {Cos[t]\ \((2 - 1.5\ Cos[2\ t]\^2)\), Sin[t]}; 
    
{a, b} = {0, 2\ π}; 


    ParametricPlot[{x[t], y[t]}, {t, a, b}, 
      PlotStyle -> {{Thickness[0.01], Blue}}, 
      AspectRatio -> Automatic, AxesLabel -> {"x", "y"}, 
      Epilog -> Text["R", {2, 3}]]; \)

Use the Gauss-Green formula to calculate
       ∫(∫_R) (x^2) y dx 
            dy.

First check to see whether the parameterization of the boundary of R is counterclockwise:

early = 
      ParametricPlot[{x[t], y[t]}, {t, a, a + 1}, 
        PlotStyle -> {{Thickness[0.01], Blue}}, 
        AspectRatio -> Automatic, AxesLabel -> {"x", "y"}];

later = 
      ParametricPlot[{x[t], y[t]}, {t, a + 1, a + 3}, 
        PlotStyle -> {{Thickness[0.02], Red}}, 
        AspectRatio -> Automatic, AxesLabel -> {"x", "y"}];

Show[early, later];

Yep; counterclockwise.
Now enter f[x, y] = (x^2) y:

\(Clear[f]; 
f[x_, y_] = x\^2\ y\)

Set n[x, y] = (∫_0^x) f[s, y] ds:

\(Clear[n, s]; 

    n[x_, y_] = ∫\_0\%x f[s, y] \[DifferentialD]s\)

Calculate
       ∫(∫_R) f[x, y] dx 
              dy = 
       ∫(∫_R) D[n[x, y], x] dx 
              dy = 
        (∫_a^b) n[x[t], y[t]] y'[t] dt:

NIntegrate[n[x[t], y[t]] y'[t], {t, a, b}]

No sweat.

B.4.d)

What else is the Gauss-Green formula good for?

If you're going on to vector calculus, you'll see that the Gauss-Green formula is at the base of a lot of the measurement procedures you'll learn.
Even if you're not going on to vector calculus, you'll see it used for a couple of snazzy things in this lesson. 
Stay tuned.

B.5)  An indication of some of the ideas behind the Gauss-Green formula

B.5.a)

What are some of the ideas behind the Gauss-Green formula?

The Gauss-Green formula says this: 
When you have a parameterization {x[t], y[t]} of the boundary of a region R in the xy-plane with the extra feature that as t advances from a to b,  the points {x[t], y[t]} advance around the boundary of a R exactly once in the counterclockwise way, then you are guaranteed that
       ∫(∫_R)(∂_x) n[x, y] - 
        (∂_y) m[x, y] dx dy
        = 
        (∫_a^b) m[x[t], y[t]] x'[t] + 
          n[x[t], y[t]] y'[t] dt.
Actually the Gauss-Green formula comes from a pair of formulas:
→   ∫(∫_R)(∂_x) n[x, y] dx 
              dy = 
        (∫_a^b) n[x[t], y[t]] y'[t] dt,
and
→   -
          ∫(∫_R)(∂_y) m[x, y] 
                dx dy = 
        (∫_a^b) m[x[t], y[t]] x'[t] dt.
You might say both of them are manifestations of the fundamental formula.
Here's the idea behind the formula
       -
          ∫(∫_R)(∂_y) m[x, y] 
                dx dy = 
        (∫_a^b) m[x[t], y[t]] x'[t] dt,
and it's going to be a hard pill to swallow:
Suppose R is the region between the curves y = high[x] and y = low[x] as follows:

c = -1; d = 5; \(Clear[high, low, x]; 

    high[x_] = 
      \((2 + 1\/2\ Sin[x\/3])\)\ \@\(1 - \((\(x - 2\)\/3)\)\^2\) + 2 + x\/8; 
    
\(low[x_] = 
      \(-\((3 + 1\/3\ Sin[1 - x\/5])\)\)\ \@\(1 - \((\(x - 2\)\/3)\)\^2\) + 
        2 + x\/8; \)\)Rplot = 
      Plot[{high[x], low[x]}, {x, c, d}, PlotRange -> All, 
        PlotStyle -> Blue, Ticks -> None, 
        AxesLabel -> {"x", "y"}, 
        DisplayFunction -> Identity]; leftright = 
      Graphics[{Line[{{c, -0.5}, {c, 3}}], Text["c", {c, 0}, {2, 1}], 
          Line[{{d, -0.5}, {d, 4}}], 
          Text["d", {d, 0}, {-2, 1}]}]; lrpoints = 
      Graphics[{PointSize[0.02], Point[{c, low[c]}], Point[{d, low[d]}]}]; 
    hllabels = 
      Graphics[{
          Text["y = high[x]", {((c + d)/2), high[((c + d)/2)]}, {0, 
              -2}], 
          Text["y = low[x]", {((c + d)/2), low[((c + d)/2)]}, {0, 
              +2}]}]; highlowsetup = 
      Show[Rplot, lrpoints, hllabels, leftright, 
        PlotRange -> {{-2, 6}, {-2.5, 5}}, 
        DisplayFunction -> $DisplayFunction];

Take any function m[x, y] and calculate
       ∫(∫_R)(∂_y) m[x, y] dx 
            dy 
as follows:
       ∫(∫_R)(∂_y) m[x, y] dx 
            dy
        = 
        ∫(∫_R)(∂_y) m[x, y] dy 
              dx
        = 
        (∫_c^d)(∫_(low[x] )^(high[x])) 
            (∂_y) m[x, y] dy 
              dx
        = 
        (∫_c^d) m[x, high[x]] - 
          m[x, low[x]] dx.
Put this aside for a little while.
A way to parameterize the boundary curve of R in accordance with the restrictions laid down in the Gauss-Green formula is to put
       {x[t], y[t]} = {t, low[t]} for c ≤ t ≤ d.
This starts on the far left of the low curve and brings you over to the far right along the low curve IN THE COUNTERCLOCKWISE fashion. Continue on the high curve from far right to far left by putting
       {x[t], y[t]} = {2  d - t, high[2  d - t]}
for d ≤ t ≤ (2  d - c).
Put a = c and b = 2  d - c and note that as t advances from a to b, {x[t], y[t]} goes around the boundary curve exactly once and it goes in the counterclockwise fashion.
Now
       (∫_a^b) m[x[t], y[t]] x'[t] dt = 
        (∫_c^(2  d - c))m[x[t], y[t]] x'[t] 
            dt
        = 
        (∫_c^d) m[x[t], y[t]] x'[t] dt + 
          (∫_d^(2  d - c))m[x[t], y[t]] x'[t] 
              dt
        = 
        (∫_c^d) m[t, low[t]] dt + 
          (∫_d^(2  d - c))m [2  d - t, 
              high[2  d - t] (-1) dt.
The substitution u = 2  d - t transforms
       (∫_d^(2  d - c))m [2  d - t, 
          high[2  d - t] (-1) dt
into
       (∫_d^c) m[u, high[u]] dt = 
        -(∫_c^d) m[u, high[u]] dt.
So
       (∫_a^b) m[x[t], y[t]] x'[t] dt
        = 
        (∫_c^d) m[t, low[t]] dt - 
          (∫_c^d) m[u, high[u]] du
        = 
        (∫_c^d) m[x, low[x]] dx - 
          (∫_c^d) m[x, high[x]] dx.
But from above,
       ∫(∫_R) D[m[x, y], y] dx 
            dy
        = 
        (∫_c^d) m[x, high[x]] - 
          m[x, low[x]] dx.
Shazam! 
The inescapable conclusion is:
       ∫(∫_R)(-D[m[x, y], y]) dx 
            dy
        = 
        (∫_a^b) m[x[t], y[t]] x'[t] dt.

To see the idea behind the companion formula
       ∫(∫_R)(D[n[x, y], x]) dx 
            dy
        = 
        (∫_a^b) n[x[t], y[t]] y'[t] dt,
do the same thing as above but modify it at the start by going with left[y] and right[y] as in the plot below; then integrate first with respect to x and second with respect to y:

c = -1; d = 4.7; highlow = 
      Graphics[{Line[{{-1, c}, {5, c}}], Text["c", {0, c}, {2, 2}], 
          Line[{{-1, d}, {5, d}}], Text["d", {0, d}, {-2, 2}]}]; 
    hlpoints = 
      Graphics[{PointSize[0.02], Point[{1.6, low[1.6]}], 
          Point[{2.77, high[2.77]}]}]; hllabels = {Graphics[Text["x = left[y]", {-0.7, 3.55}]], 
        Graphics[Text["x = right[y]", {4.6, 2.5}]]}; highlowsetup = 
      Show[Rplot, hlpoints, hllabels, highlow, 
        PlotRange -> {{-2, 6}, {-2.5, 5}}, 
        DisplayFunction -> $DisplayFunction];

At this point, you have the idea behind the formulas
       ∫(∫_R)(-(∂_y) m[x, y]) 
            dx dy
        = 
        (∫_a^b) m[x[t], y[t]] x'[t] dt
and
       ∫(∫_R)(∂_x) n[x, y] dx 
            dy
        = 
        (∫_a^b) n[x[t], y[t]] y'[t] dt.
Add both formulas together to get the full-blown Gauss-Green formula
       ∫(∫_R)(∂_x) n[x, y] - 
        (∂_y) m[x, y] dx dy
        = 
        (∫_a^b) m[x[t], y[t]] x'[t] + 
          n[x[t], y[t]] y'[t] dt.
'Nuff said.