Lecture March 1 -- Approximation of integrals by Riemann sums

Initialization cell.  Execute this group of cells before beginning

Left endpoint approximation

\(\(\(leftrect[n_] := 
        Block[{X, dX, rectangles, boundaries, RightRect, \ leftsum, \ label, 
            aloc, \ bloc, floc}, 
aloc\  = \ 
            a; 
bloc\  = \ b; 
floc\  = \ 
            f; 
functionplot := 
            Plot[floc[x], {x, a, b}, DisplayFunction -> Identity, 
              PlotRange -> All]; 
RightRect[g_, x_, 
              d_] := Line[{{x, 0}, {x + d, 0}, {x + d, g[x]}, {x, g[x]}, {x, 
                  0}}]; 
dX = \(bloc - aloc\)\/n; 
          X = Table[a + \(i\ \((bloc - aloc)\)\)\/n, {i, 0, n}]; 
          rectangles := 
            Graphics[{GrayLevel[0.9], 
                Table[Rectangle[{X[[i]], 
                      0}, {X[[i + 1]], 
                      f[X[[i]]]}], {i, 
                    1, n}]}]; 
          boundaries := 
            Graphics[
              Table[RightRect[floc, 
                  X[[j]], dX], {j, 1, 
                  n}]]; leftsum = 
            N[∑\+\(i = 1\)\%n floc[
                    X[[i]]]\ dX, 10]; 
          label = leftsum\ "=  sum of  rectangle areas "; 
          Show[{functionplot, \ rectangles, boundaries, \ functionplot}, 
            PlotLabel -> label, 
            Frame\  -> \ 
              True, \ \ DisplayFunction -> $DisplayFunction];];\)\(
\)
    \)\)
\(\(leftrect[f_, \ a_, \ b_, \ n_]\  := \ 
        Module[{X, dX, rectangles, boundaries, RightRect, \ leftsum, \ label, 
            floc, \ aloc, \ bloc}, 
aloc\  = \ 
            a; 
bloc\  = \ b; 
floc\  = \ 
            f; 
functionplot := 
            Plot[floc[x], {x, a, b}, DisplayFunction -> Identity, 
              PlotRange -> All]; 
RightRect[g_, x_, 
              d_] := Line[{{x, 0}, {x + d, 0}, {x + d, g[x]}, {x, g[x]}, {x, 
                  0}}]; 
dX = \(bloc - aloc\)\/n; 
          X = Table[a + \(i\ \((bloc - aloc)\)\)\/n, {i, 0, n}]; 
          rectangles := 
            Graphics[{GrayLevel[0.9], 
                Table[Rectangle[{X[[i]], 
                      0}, {X[[i + 1]], 
                      f[X[[i]]]}], {i, 
                    1, n}]}]; 
          boundaries := 
            Graphics[
              Table[RightRect[floc, 
                  X[[j]], dX], {j, 1, 
                  n}]]; leftsum = 
            N[∑\+\(i = 1\)\%n floc[
                    X[[i]]]\ dX, 10]; 
          label = leftsum\ "=  sum of  rectangle areas "; 
          Show[{functionplot, \ rectangles, boundaries, \ functionplot}, 
            PlotLabel -> label, 
            Frame\  -> \ 
              True, \ \ DisplayFunction -> $DisplayFunction];];\)


\)
leftsum[n_]  :=  
        Block[{X, dX, rectangles, boundaries, 
            RightRect}, 
aloc  =  
            a; 
bloc  =  
            b; 
dX = ((bloc - aloc)/n); 
X = 
            Table[a + ((i (bloc - aloc))/n), {i, 0, 
                n}]; 

N[(∑_(i = 1)^n) f[X[[i]]] dX, 
            10]];

    
leftsum[f_,  a,  b_,  n_]  := 
        Block[{X, dX, rectangles, boundaries, 
            RightRect}, 
aloc  =  
            a; 
bloc  =  
            b; 
dX = ((bloc - aloc)/n); 
X = 
            Table[a + ((i (bloc - aloc))/n), {i, 0, 
                n}]; 

N[(∑_(i = 1)^n) f[X[[i]]] dX, 
            10]];

Midpoint approximation

Clear[midrect,  midsum];

    
\(\(\(midrect[n_] := 
        Block[{X, dX, rectangles, boundaries, MidRect, \ midsum, \ label, \ 
            aloc, \ bloc, \ floc}, 
aloc\  = \ 
            a; 
bloc\  = \ b; 
floc\  = \ 
            f; 
functionplot := 
            Plot[floc[x], {x, a, b}, PlotRange\  -> \ All, \ 
              DisplayFunction -> Identity]; 
          MidRect[g_, x_, d_] := 
            Line[{{x, 0}, {x + d, 0}, {x + d, g[x + d\/2]}, {x, 
                  g[x + d\/2]}, {x, 0}}]; dX = \(bloc - aloc\)\/n; 
          X = Table[aloc + \(i\ \((bloc - aloc)\)\)\/n, {i, 0, n}]; 
          rectangles := 
            Graphics[{GrayLevel[0.9], 
                Table[Rectangle[{X[[i]], 
                      0}, {X[[i + 1]], 
                      floc[X[[i]] + 
                          dX\/2]}], {i, 1, n}]}]; 
          boundaries := 
            Graphics[
              Table[MidRect[floc, 
                  X[[j]], dX], {j, 1, 
                  n}]]; midsum = 
            N[∑\+\(i = 1\)\%n floc[
                    X[[i]] + dX\/2]\ dX, 
              10]; label = midsum\ "= sum  of rectangle areas"; 
          Show[{functionplot\ , \ rectangles, boundaries, functionplot}, 
            PlotLabel -> label, 
\ 
            DisplayFunction -> $DisplayFunction];];\)\(
\)
    \)\)
midrect[f_,  a_,  b_,  n_] := 
        Block[{X, dX, rectangles, boundaries, MidRect,  midsum,  label,  
            aloc,  bloc,  floc}, 
aloc  =  
            a; 
bloc  =  b; 
floc  =  
            f; 
functionplot := 
            Plot[floc[x], {x, a, b}, PlotRange  ->  All,  
              DisplayFunction -> Identity]; 
          MidRect[g_, x_, d_] := 
            Line[{{x, 0}, {x + d, 0}, {x + d, g[x + (d/2)]}, {x, 
                  g[x + (d/2)]}, {x, 0}}]; dX = ((bloc - aloc)/n); 
          X = Table[aloc + ((i (bloc - aloc))/n), {i, 0, n}]; 
          rectangles := 
            Graphics[{GrayLevel[0.9], 
                Table[Rectangle[{X[[i]], 
                      0}, {X[[i + 1]], 
                      floc[X[[i]] + 
                          (dX/2)]}], {i, 1, n}]}]; 
          boundaries := 
            Graphics[
              Table[MidRect[floc, 
                  X[[j]], dX], {j, 1, 
                  n}]]; midsum = 
            N[(∑_(i = 1)^n) floc[
                    X[[i]] + (dX/2)] dX, 
              10]; label = midsum "= sum  of rectangle areas"; 
          Show[{functionplot ,  rectangles, boundaries, functionplot}, 
            PlotLabel -> label, 
 
            DisplayFunction -> $DisplayFunction];];
midsum[n_]  :=  
        Block[{X, dX, rectangles, boundaries, 
            RightRect}, 
aloc  =  
            a; 
bloc  =  
            b; 
dX = ((bloc - aloc)/n); 
X = 
            Table[a + ((i (bloc - aloc))/n), {i, 0, 
                n}]; 

N[(∑_(i = 1)^n) f[X[[i]]  +  dX/2] dX, 
            10]];

    
midsum[f_,  a_,  b_,  n_]  := 
        Block[{X, dX, rectangles, boundaries, 
            RightRect}, 
aloc  =  
            a; 
bloc  =  
            b; 
dX = ((bloc - aloc)/n); 
X = 
            Table[a + ((i (bloc - aloc))/n), {i, 0, 
                n}]; 

N[(∑_(i = 1)^n) f[X[[i]] +  dX/2] dX, 
            10]];

Right endpoint approximation

\(\(rightrect[n_] := 
        Block[{X, dX, rectangles, boundaries, RightRect, \ floc, \ aloc, \ 
            bloc, h}, 
floc\  = \ 
            f; 
aloc\  = \ a; 
bloc\  = \ 
            b; 
h\  = \ 
            b - a; 
functionplot := 
            Plot[floc[x], {x, a, b}, DisplayFunction -> Identity, 
              PlotRange -> All]; 
          RightRect[g_, x_, d_] := 
            Line[{{x, 0}, {x + d, 0}, {x + d, g[x + d]}, {x, g[x + d]}, {x, 
                  0}}]; dX = h\/n; 
          X = Table[a + \(i\ \((h)\)\)\/n, {i, 0, n}]; 
          rectangles := 
            Graphics[{GrayLevel[0.9], 
                Table[Rectangle[{X[[i]], 
                      0}, {X[[i + 1]], 
                      floc[X[[i + 
                            1]]]}], {i, 1, n}]}]; 
          boundaries := 
            Graphics[
              Table[RightRect[floc, 
                  X[[j]], dX], {j, 1, 
                  n}]]; rightsum = 
            N[∑\+\(i = 1\)\%n floc[
                    X[[i + 1]]]\ dX, 
              10]; label = rightsum\ "=  sum of area of rectangles "; 
          Show[{functionplot, \ rectangles, boundaries, functionplot}, 
            PlotLabel -> 
              label, 
DisplayFunction -> $DisplayFunction];];\)\)
\(\(rightrect[f_, \ a_, \ b_, \ n_] := 
        Block[{X, dX, rectangles, boundaries, RightRect, \ floc, \ aloc, \ 
            bloc, h}, 
floc\  = \ 
            f; 
aloc\  = \ a; 
bloc\  = \ 
            b; 
h\  = \ 
            b - a; 
functionplot := 
            Plot[floc[x], {x, a, b}, DisplayFunction -> Identity, 
              PlotRange -> All]; 
          RightRect[g_, x_, d_] := 
            Line[{{x, 0}, {x + d, 0}, {x + d, g[x + d]}, {x, g[x + d]}, {x, 
                  0}}]; dX = h\/n; 
          X = Table[a + \(i\ \((h)\)\)\/n, {i, 0, n}]; 
          rectangles := 
            Graphics[{GrayLevel[0.9], 
                Table[Rectangle[{X[[i]], 
                      0}, {X[[i + 1]], 
                      floc[X[[i + 
                            1]]]}], {i, 1, n}]}]; 
          boundaries := 
            Graphics[
              Table[RightRect[floc, 
                  X[[j]], dX], {j, 1, 
                  n}]]; rightsum = 
            N[∑\+\(i = 1\)\%n floc[
                    X[[i + 1]]]\ dX, 
              10]; label = rightsum\ "=  sum of area of rectangles "; 
          Show[{functionplot, \ rectangles, boundaries, functionplot}, 
            PlotLabel -> 
              label, 
DisplayFunction -> $DisplayFunction];];\)\)

Trapezoid approximation

trapezoids[n_] := 
        Block[{X, dX, Trap, regions, boundaries}, 
          functionplot := 
            Plot[f[x], {x, a, b}, DisplayFunction -> Identity, 
              PlotRange -> All]; dX = ((b - a)/n); 
          X = Table[a + ((i (b - a))/n), {i, 0, n}]; 
          Trap[g_, x_, d_] := 
            Line[{{x, 0}, {x + d, 0}, {x + d, g[x + d]}, {x, g[x]}, {x, 0}}]; 
          regions := 
            Graphics[{GrayLevel[0.9], 
                Table[Polygon[{{X[[i]], 
                        0}, {X[[i]], 
                        f[X[[
                            i]]]}, {X[[
                          i + 1]], 
                        f[X[[i + 
                              1]]]}, {X[[
                          i + 1]], 0}}], {i, 1, n}]}]; 
          boundaries := 
            Graphics[
              Table[Trap[f, X[[j]], 
                  dX], {j, 1, n}]]; 
          trapsum = 
            N[(1/2) dX (f[a] + 
                    2 (∑_(i = 2)^n) f[
                          X[[i]]] + 
                    f[b]), 10]; 
          label = trapsum "= sum of area of trapezoids"; 
          Show[regions, boundaries, functionplot, PlotLabel -> label, 
            DisplayFunction -> $DisplayFunction];];

Numerical rules

LeftSum[n_] := 
      Block[{dX}, dX = ((b - a)/n); 
        N[(∑_(i = 0)^(n - 1))dX f[a + i dX], 10]]
MidpointSum[n_] := 
      Block[{dX}, dX = ((b - a)/n); 
        N[(∑_(i = 0)^(n - 1))dX f[a + (i + 0.5) dX], 
          10]]
RightSum[n_] := 
      Block[{dX}, dX = ((b - a)/n); 
        N[(∑_(i = 1)^n) dX f[a + i dX], 10]]
TrapezoidSum[n_] := 
      Block[{dX}, dX = ((b - a)/n); 
        N[(∑_(i = 0)^(n - 1))dX 0.5 (f[a + i dX] + 
                  f[a + (i + 1) dX]), 10]]

2d integration stuff

Needs["Graphics`Colors`"];
xFace[f_, x_, c_, d_] := 
        Module[{s, yy, zz, face, intersection}, 
zz = 
            s f[x, yy]; 
face = 
            ParametricPlot3D[{x, yy, zz}, {s, 0, 1}, {yy, c, d}, 
              PlotPoints -> {2, 30}, Boxed -> False, 
              Axes -> False, 
              DisplayFunction -> Identity]; 
face = 
            Graphics3D[{EdgeForm[], First[face]}]; 
zz = 
            f[x, yy]; 
intersection = 
            ParametricPlot3D[{x, yy, zz}, {yy, c, d}, Boxed -> False, 
              Axes -> False, 
              DisplayFunction -> Identity]; 
Return[
            Show[{face, intersection, 
                Graphics3D[{Line[{{x, d, f[x, d]}, {x, d, 0}, {x, c, 0}, {x, 
                          c, f[x, c]}}]}]}, 
              DisplayFunction -> 
                Identity]];
];

    
yFace[f_, a_, b_, y_] := 
        Module[{s, xx, zz, face, intersection}, 
zz = 
            s f[xx, y]; 
face = 
            ParametricPlot3D[{xx, y, zz}, {s, 0, 1}, {xx, a, b}, 
              PlotPoints -> {2, 30}, Boxed -> False, 
              Axes -> False, 
              DisplayFunction -> Identity]; 
face = 
            Graphics3D[{EdgeForm[], First[face]}]; 
zz = 
            f[xx, y]; 
intersection = 
            ParametricPlot3D[{xx, y, zz}, {xx, a, b}, Boxed -> False, 
              Axes -> False, 
              DisplayFunction -> Identity]; 
Return[
            Show[{face, intersection, 
                Graphics3D[{Line[{{b, y, f[b, y]}, {b, y, 0}, {a, y, 0}, {a, 
                          y, f[a, y]}}]}]}, 
              DisplayFunction -> 
                Identity]];
];

    
fDxDyColumn[f_, x_, y_, dx_] := 
        Module[{z, column}, 
z = 
            f[x  +  dx/2, y  +  dx/2]; 
Return[
            Graphics3D[{Blue, 
                Line[{{x, y, z}, {x + dx, y, z}, {x + dx, y + dx, z}, {x, 
                      y + dx, z}, {x, y, z}, {x + dx, y, z}, {x + dx, y, 
                      0}, {x, y, 0}, {x, y, z}, {x, y + dx, z}, {x, y + dx, 
                      0}, {x, y, 0}, {x, y, z}, {x, y + dx, z}, {x + dx, 
                      y + dx, z}, {x + dx, y + dx, 0}, {x, y + dx, 0}, {x, 
                      y + dx, z}, {x, y + dx, 0}, {x + dx, y + dx, 
                      0}, {x + dx, y, 0}, {x, y, 0}, {x, y + dx, 0}, {x + dx, 
                      y + dx, 0}, {x + dx, y, 0}, {x + dx, y, z}, {x + dx, 
                      y + dx, z}, {x, y + dx, z}, {x, y, 
                      z}
}], 
                Polygon[{{x, y, z}, {x + dx, y, z}, {x + dx, y + dx, z}, {x, 
                      y + dx, z}, {x, y, z}, {x + dx, y, z}, {x + dx, y, 
                      0}, {x, y, 0}, {x, y, z}, {x, y + dx, z}, {x, y + dx, 
                      0}, {x, y, 0}, {x, y, z}, {x, y + dx, z}, {x + dx, 
                      y + dx, z}, {x + dx, y + dx, 0}, {x, y + dx, 0}, {x, 
                      y + dx, z}, {x, y + dx, 0}, {x + dx, y + dx, 
                      0}, {x + dx, y, 0}, {x, y, 0}, {x, y + dx, 0}, {x + dx, 
                      y + dx, 0}, {x + dx, y, 0}, {x + dx, y, z}, {x + dx, 
                      y + dx, z}, {x, y + dx, z}, {x, y, 
                      z}
}]}]];
;];

    
fDySlab[f_, a_, b_, y_, dy_] := 
        Module[{xx, yy, ff, ttop}, 
ff[xx_, yy_] = 
            f[xx, y]; 

ttop = 
            Plot3D[ff[xx, yy], {xx, a, b}, {yy, y, y + dy}, 
              Ticks -> None, PlotPoints -> {30, 2}, 
              Mesh -> False, 
              DisplayFunction -> 
                Identity]; 

Return[
            Show[{ttop, 
                Graphics3D[{Polygon[{{a, y, f[a, y]}, {a, y, 0}, {a, y + dy, 
                          0}, {a, y + dy, f[a, y]}}], 
                    Line[{{a, y, f[a, y]}, {a, y, 0}, {a, y + dy, 0}, {a, 
                          y + dy, f[a, y]}}], 
                    Polygon[{{b, y, f[b, y]}, {b, y, 0}, {b, y + dy, 0}, {b, 
                          y + dy, f[b, y]}}], 
                    Line[{{b, y, f[b, y]}, {b, y, 0}, {b, y + dy, 0}, {b, 
                          y + dy, f[b, y]}}]}], yFace[ff, a, b, y], 
                yFace[ff, a, b, y + dy]}, 
              DisplayFunction -> Identity]]];

    
fDxSlab[f_, x_, dx_, c_, d_] := 
        Module[{xx, yy, ff, ttop}, 
ff[xx_, yy_] = 
            f[x, yy]; 

ttop = 
            Plot3D[ff[xx, yy], {xx, x, x + dx}, {yy, c, d}, 
              Ticks -> None, PlotPoints -> {2, 30}, 
              Mesh -> False, 
              DisplayFunction -> 
                Identity]; 

Return[
            Show[{ttop, xFace[ff, x, c, d], xFace[ff, x + dx, c, d], 
                Graphics3D[{Polygon[{{x, c, f[x, c]}, {x, c, 0}, {x + dx, c, 
                          0}, {x + dx, c, f[x, c]}}], 
                    Line[{{x, c, f[x, c]}, {x, c, 0}, {x + dx, c, 
                          0}, {x + dx, c, f[x, c]}}], 
                    Polygon[{{x, d, f[x, d]}, {x, d, 0}, {x + dx, d, 
                          0}, {x + dx, d, f[x, d]}}], 
                    Line[{{x, d, f[x, d]}, {x, d, 0}, {x + dx, d, 
                          0}, {x + dx, d, f[x, d]}}]}]}, 
              DisplayFunction -> 
                Identity]]];

Illustration of limits of Riemann sums approximating an integral, for a function with both positive and negative values.  (The integral represents a sum of signed areas.)

f[x_]  :=  
      x^3  -  4  x^2  + 5

    Plot[f[x],  {x,  -3,  7}]

midrect[f,  -2,  4,  20]
midrect[f,  -2,  4,  100]
midrect[f,  -2,  4,  500]

Integrate[f[x],  {x, -2, 4}]

Riemann sum approximation for an integral with respect to area.  (The integral represents a volume.)

f[x_, y_] := 4 - (x^2) - ((y^2)/3);
z = f[x, y];
a = -1; b = (3/2);
c = 0; d = 2;
\(\(\(viewArrow = {5, \(-3\), 2};\)\(
\)
    \)\)
dx = 1/2;

top = 
        Plot3D[z, {x, a, b}, {y, c, d}, BoxRatios -> Automatic, 
          ViewPoint -> viewArrow, Ticks -> None, 
          AxesLabel -> {"x", "y", "z"}, Mesh -> False, 
          DisplayFunction -> Identity];
faces = 
        Show[{xFace[f, a, c, d], xFace[f, b, c, d], yFace[f, a, b, c], 
            yFace[f, a, b, d]}, 
          PlotRange -> 1.01 {{a, b}, {c, d}, {0, 4}}, 
          ViewPoint -> viewArrow, 
          DisplayFunction -> Identity];
vol = 
      Show[{faces, top}, PlotRange -> 1.01 {{a, b}, {c, d}, {0, 4}}, 
        ViewPoint -> viewArrow, DisplayFunction -> Identity]; 
    approx = Table[
        fDxDyColumn[f, x, y, dx], {x, a, b - dx, dx}, {y, c, d - dx, 
          dx}];
approx = 
        Show[approx, ViewPoint -> viewArrow, 
          DisplayFunction -> Identity];
Show[GraphicsArray[{vol, approx}], 
        DisplayFunction -> $DisplayFunction];

dx = 1/6;

top = 
        Plot3D[z, {x, a, b}, {y, c, d}, BoxRatios -> Automatic, 
          ViewPoint -> viewArrow, Ticks -> None, 
          AxesLabel -> {"x", "y", "z"}, Mesh -> False, 
          DisplayFunction -> Identity];
faces = 
        Show[{xFace[f, a, c, d], xFace[f, b, c, d], yFace[f, a, b, c], 
            yFace[f, a, b, d]}, 
          PlotRange -> 1.01 {{a, b}, {c, d}, {0, 4}}, 
          ViewPoint -> viewArrow, 
          DisplayFunction -> Identity];
vol = 
      Show[{faces, top}, PlotRange -> 1.01 {{a, b}, {c, d}, {0, 4}}, 
        ViewPoint -> viewArrow, DisplayFunction -> Identity]; 
    approx = Table[
        fDxDyColumn[f, x, y, dx], {x, a, b - dx, dx}, {y, c, d - dx, 
          dx}];
approx = 
        Show[approx, ViewPoint -> viewArrow, 
          DisplayFunction -> Identity];
Show[GraphicsArray[{vol, approx}], 
        DisplayFunction -> $DisplayFunction];