(************** Content-type: application/mathematica **************

                    Mathematica-Compatible Notebook

This notebook can be used with any Mathematica-compatible
application, such as Mathematica, MathReader or Publicon. The data
for the notebook starts with the line containing stars above.

To get the notebook into a Mathematica-compatible application, do
one of the following:

* Save the data starting with the line of stars above into a file
  with a name ending in .nb, then open the file inside the
  application;

* Copy the data starting with the line of stars above to the
  clipboard, then use the Paste menu command inside the application.

Data for notebooks contains only printable 7-bit ASCII and can be
sent directly in email or through ftp in text mode.  Newlines can be
CR, LF or CRLF (Unix, Macintosh or MS-DOS style).

NOTE: If you modify the data for this notebook not in a Mathematica-
compatible application, you must delete the line below containing
the word CacheID, otherwise Mathematica-compatible applications may
try to use invalid cache data.

For more information on notebooks and Mathematica-compatible 
applications, contact Wolfram Research:
  web: http://www.wolfram.com
  email: info@wolfram.com
  phone: +1-217-398-0700 (U.S.)

Notebook reader applications are available free of charge from 
Wolfram Research.
*******************************************************************)

(*CacheID: 232*)


(*NotebookFileLineBreakTest
NotebookFileLineBreakTest*)
(*NotebookOptionsPosition[     18304,        480]*)
(*NotebookOutlinePosition[     19043,        505]*)
(*  CellTagsIndexPosition[     18999,        501]*)
(*WindowFrame->Normal*)



Notebook[{
Cell["\<\
Go through the entire lab and turn in the two You Try It portions, together \
with a 1-2 page summary of the substance of the lab. 
Be sure to contrast the error patterns you observed when you sampled more \
rather than fewer points and repeated your sampling many times. Also contrast \
your error patterns with those of your classmates.\
\>", "Text"],

Cell[CellGroupData[{

Cell["\<\
Take Your Chances: Try the Monte Carlo Technique for Numerical Integration in \
Three Dimensions  \
\>", "Title"],

Cell[TextData[StyleBox["Chapter 12, Section 1",
  FontFamily->"Arial",
  FontSize->16,
  FontWeight->"Bold"]], "Text"],

Cell[CellGroupData[{

Cell["Introduction", "Section"],

Cell["\<\
OBJECTIVE: Learn the Monte Carlo method for approximating a multiple integral \
that cannot be integrated symbolically.\
\>", "Text"],

Cell["\<\
How can you use a game of chance to evaluate a multiple integral? That is \
just what you will do with this project. You will generate random points \
within a fixed region and then estimate the volume of the desired portion by \
considering the percentage of random points that fall within the boundaries \
of the desired portion. Since this will only estimate the exact volume, you \
will also explore the accuracy of this method.\
\>", "Text"],

Cell[CellGroupData[{

Cell["Technology Guidelines", "Subsection",
  CellDingbat->"\[LightBulb]"],

Cell[TextData[{
  StyleBox["NOTE:  If you have just finished a module, restart ",
    CellFrame->True,
    Background->None],
  StyleBox["Mathematica",
    CellFrame->True,
    FontSlant->"Italic",
    Background->None],
  StyleBox[
  " before executing a new module.\nTO OPEN CELLS, put your cursor on the \
right cell bracket and double click.",
    CellFrame->True,
    Background->None],
  "\nTO STOP AN EXECUTION\n\tSelect the ",
  StyleBox["Kernel",
    FontSlant->"Italic"],
  " pull-down menu and click on ",
  StyleBox["Abort Evaluation.\n",
    FontSlant->"Italic"],
  "ORDER OF EXECUTION\n\tExecute cells in the order given. Do not skip any \
Input cells within a given notebook.\nSAVING NOTEBOOKS\n\tYou can save \
anytime to any directory you choose, and it is wise to save often.\n\t\
However, before you do your final save, delete all your output by selecting \
the \n\t ",
  StyleBox["Delete All Output",
    FontSlant->"Italic"],
  " selection under the ",
  StyleBox["Kernel",
    FontSlant->"Italic"],
  " pull-down menu.\nEXPERIENCING MAJOR PROBLEMS\n\tSave if appropriate, and \
then shut down ",
  StyleBox["Mathematica",
    FontSlant->"Italic"],
  " and start it up again."
}], "Text"]
}, Closed]]
}, Closed]],

Cell[CellGroupData[{

Cell["Part I: The Monte Carlo Method", "Section"],

Cell[TextData[{
  "The Monte Carlo method is an example of a numerical integration technique. \
In this method, the volume under a surface is completely enclosed in a \
rectangular box and points within that box are randomly selected. The number \
of points that land under the surface is determined and the volume under the \
surface to the bottom of the box is approximated as the percentage of the \
generated points lying under the surface times the volume of the entire box. \
In the example that follows, 1000 points are generated randomly to estimate \
the volume. Note that the sample function used here cannot be integrated \
symbolically, so we are using ",
  StyleBox["Mathematica",
    FontSlant->"Italic"],
  "'s numerical integrator, which uses a different approximation technique, \
to come up with an answer that we would expect to be close to the actual one.\
\n\nIn utilizing the ",
  StyleBox["Random",
    FontWeight->"Bold"],
  " command in ",
  StyleBox["Mathematica",
    FontSlant->"Italic"],
  ", the default ",
  StyleBox["Random[]",
    FontWeight->"Bold"],
  " selects numbers between 0 and 1; other selections are specified. The \
following commands will make the dimesions of the box selected 2 by 3 by 1, \
giving a volume of 6. When applying a function of your own to this procedure, \
be certain to carefully construct the box in which you enclose your desired \
volume. In choosing a function, be careful to choose one that is not negative \
in the specifed domain, and note that you are finding the volume between the \
surface and the ",
  StyleBox["x-y",
    FontSlant->"Italic"],
  " plane. We begin by defining the function, specifying the region over \
which we will generate points, and then looking at the graph of the function \
and the region"
}], "Text"],

Cell[BoxData[{
    \(Off[General::spell]\ \), "\n", 
    \(Off[General::spell1]\ \), "\n", 
    \(Clear[x, y, z, f]\), "\n", 
    \(f[x_, y_] := E\^\(\(-2\) x\^2 - y\^2\)\), "\n", 
    \(\(xmin = 0;\)\), "\n", 
    \(\(xmax = 1.5;\)\), "\n", 
    \(\(ymin = 0;\)\), "\n", 
    \(\(ymax = 2;\)\), "\n", 
    \(\(zmin = 0;\)\), "\n", 
    \(\(zmax = 1;\)\), "\n", 
    \(\(pf = 
        Plot3D[f[x, y], {x, xmin, xmax}, {y, ymin, ymax}, 
          PlotRange -> All];\)\), "\n", 
    \(\(actual = 
        NIntegrate[f[x, y], {x, xmin, xmax}, {y, ymin, ymax}];\)\), "\n", 
    \(Print["\<actual volume = \>", actual]\)}], "Input"],

Cell[TextData[{
  "Now we will generate random points in our specified region and check to \
see whether or not they are under our surface. The points in ",
  StyleBox["listin",
    FontWeight->"Bold"],
  " are under the surface in the volume we are computing and they will be \
plotted in red. The points in ",
  StyleBox["listout",
    FontWeight->"Bold"],
  " are not under the surface and they will be plotted in black. The scatter \
plot command requires the following package which must be read in first."
}], "Text"],

Cell[BoxData[
    \(<< Graphics`Graphics3D`\)], "Input"],

Cell[BoxData[{
    \(\(numberofpoints = 1000;\)\), "\n", 
    \(\(volumeofbox = \((xmax - xmin)\)\ \((ymax - ymin)\)\ \((zmax - 
              zmin)\);\)\), "\n", 
    \(\(count = 0;\)\), "\n", 
    \(\(listout = {};\)\), "\n", 
    \(\(listin = {};\)\), "\n", 
    \(Do[{Clear[x, y, z], x = Random[Real, {xmin, xmax}], 
        y = Random[Real, {ymin, ymax}], 
        z = Random[Real, {zmin, zmax}], \n\t\t\ \ \ \ \ \ If[
          z < f[x, y], {count = count + 1, AppendTo[listin, {x, y, z}]}, 
          AppendTo[listout, {x, y, z}]]}, {numberofpoints}]\), "\n", 
    \(\(pout = 
        ScatterPlot3D[listout, DisplayFunction -> Identity];\)\), "\n", 
    \(\(pin = 
        ScatterPlot3D[listin, PlotStyle -> RGBColor[1, 0, 0], 
          DisplayFunction -> Identity];\)\), "\n", 
    \(\(Show[pin, pout, pf, DisplayFunction -> $DisplayFunction];\)\), "\n", 
    \(Print["\<actual = \>", actual]\), "\n", 
    \(\(estimate = N[volumeofbox\ \ count/numberofpoints];\)\), "\n", 
    \(Print["\<estimate = \>", estimate]\), "\n", 
    \(\(error = estimate - actual;\)\), "\n", 
    \(Print["\<error = \>", error]\), "\n", 
    \(\(percenterror = error/actual;\)\), "\n", 
    \(Print["\<percenterror = \>", percenterror]\)}], "Input"],

Cell["\<\
How does your error compare to that of one of your classmates or to a second \
calculation you produce? Note what happens when you increase the number of \
points you use for the calculation.\
\>", "Text"],

Cell[BoxData[{
    \(Clear[x, y, z]\), "\n", 
    \(\(numberofpoints = 3000;\)\), "\n", 
    \(\(count = 0;\)\), "\n", 
    \(\(listout = {};\)\), "\n", 
    \(\(listin = {};\)\), "\n", 
    \(Do[{Clear[x, y, z], x = Random[Real, {xmin, xmax}], 
        y = Random[Real, {ymin, ymax}], 
        z = Random[Real, {zmin, zmax}], \n\t\t\ \ \ \ \ \ If[
          z < f[x, y], {count = count + 1, AppendTo[listin, {x, y, z}]}, 
          AppendTo[listout, {x, y, z}]]}, {numberofpoints}]\), "\n", 
    \(\(pout = 
        ScatterPlot3D[listout, DisplayFunction -> Identity];\)\), "\n", 
    \(\(pin = 
        ScatterPlot3D[listin, PlotStyle -> RGBColor[1, 0, 0], 
          DisplayFunction -> Identity];\)\), "\n", 
    \(\(Show[pin, pout, pf, DisplayFunction -> $DisplayFunction];\)\), "\n", 
    \(Print["\<actual = \>", actual]\), "\n", 
    \(\(estimate = N[volumeofbox\ \ count/numberofpoints];\)\), "\n", 
    \(Print["\<estimate = \>", estimate]\), "\n", 
    \(\(error = estimate - actual;\)\), "\n", 
    \(Print["\<error = \>", error]\), "\n", 
    \(\(percenterror = error/actual;\)\), "\n", 
    \(Print["\<percenterror = \>", percenterror]\)}], "Input"],

Cell["\<\
Does this method seem very accurate?  We will analyze the error involved in \
such an approximation.\
\>", "Text"]
}, Closed]],

Cell[CellGroupData[{

Cell["You Try It: Part I", "Section"],

Cell[TextData[{
  "Choose your own function; just remember to make it nonnegative over the \
region of integration. The only commands you need to enter in the following \
template",
  Cell[BoxData[
      \(TraditionalForm\`-\)]],
  "your function and minimum and maximum values for each of your variables \
(leave your ",
  StyleBox["zmin",
    FontWeight->"Bold"],
  " at 0)",
  Cell[BoxData[
      \(TraditionalForm\`-\)]],
  "are in red . Then, simply execute the rest of the notebook, excluding the \
first cell that defined the function used above. \n\nYou should realize that \
the bigger your rectangular region, the more points you will need to generate \
to get an answer that is reasonably accurate. Why?\n\nWe will begin by \
plotting the function and letting ",
  StyleBox["Mathematica",
    FontSlant->"Italic"],
  " estimate the volume between the surface and the ",
  StyleBox["xy",
    FontSlant->"Italic"],
  "-plane over the domain specified. Be sure to stay above the ",
  StyleBox["xy",
    FontSlant->"Italic"],
  "-plane."
}], "Text"],

Cell[BoxData[{\(Clear[x, y, z, g]\), 
    RowBox[{
      RowBox[{\(g[x_, y_]\), "=", 
        StyleBox[\(Sin[x - y]\^2\),
          FontColor->RGBColor[1, 0, 0]]}], ";", "\n", 
      RowBox[{"xmin", "=", 
        StyleBox[\(-4\),
          FontColor->RGBColor[1, 0, 0]]}], ";", "\n", 
      RowBox[{"xmax", "=", 
        StyleBox["4",
          FontColor->RGBColor[1, 0, 0]]}], ";", "\n", 
      RowBox[{"ymin", "=", 
        StyleBox["0",
          FontColor->RGBColor[1, 0, 0]]}], ";", "\n", 
      RowBox[{"ymax", "=", 
        StyleBox["3",
          FontColor->RGBColor[1, 0, 0]]}], ";", "\n", 
      RowBox[{"zmin", "=", 
        StyleBox["0",
          FontColor->RGBColor[1, 0, 0]]}], ";", "\n", 
      RowBox[{"zmax", "=", 
        StyleBox["1",
          FontColor->RGBColor[1, 0, 0]]}], ";", "\n", 
      \(pf = Plot3D[g[x, y], {x, xmin, xmax}, {y, ymin, ymax}, 
          PlotRange -> All, PlotPoints \[Rule] 40]\), ";", "\n", 
      \(actual = NIntegrate[g[x, y], {x, xmin, xmax}, {y, ymin, ymax}]\), 
      ";", "\n", \(Print["\<actual volume = \>", actual]\)}]}], "Input"],

Cell[BoxData[{
    \(numberofpoints = 2000; \n
    volumeofbox = \((xmax - xmin)\)\ \((ymax - ymin)\)\ \((zmax - zmin)\); \n
    count = 0; \nlistout = {}; \nlistin = {}; \n
    Do[{Clear[x, y, z], x = Random[Real, {xmin, xmax}], 
        y = Random[Real, {ymin, ymax}], z = Random[Real, {zmin, zmax}], \n
        \t\t\ \ \ \ \ \ 
        If[z < g[x, y], {count = count + 1, AppendTo[listin, {x, y, z}]}, 
          AppendTo[listout, {x, y, z}]]}, {numberofpoints}]\), 
    \(pout = ScatterPlot3D[listout, DisplayFunction -> Identity]; \n
    pin = ScatterPlot3D[listin, PlotStyle -> RGBColor[1, 0, 0], 
        DisplayFunction -> Identity]; \n
    Show[pin, pout, pf, DisplayFunction -> $DisplayFunction]; \n
    Print["\<actual = \>", actual]\), 
    \(estimate = N[volumeofbox\ \ count/numberofpoints]; \n
    Print["\<estimate = \>", estimate]\), 
    \(error = estimate - actual; \nPrint["\<error = \>", error]\), 
    \(percenterror = error/actual; \n
    Print["\<percenterror = \>", percenterror]\)}], "Input"]
}, Closed]],

Cell[CellGroupData[{

Cell["Part II: Error Analysis", "Section"],

Cell["\<\
In this part, we follow the procedure above with 500 points; then, we repeat \
that procedure 200 times and look at the pattern of the errors. Because of \
the many steps involved, this computation may take several minutes.\
\>", "Text"],

Cell[BoxData[{
    \(Clear[x, y, z]\), "\n", 
    \(\(numberofpoints = 500;\)\), "\n", 
    \(\(errorlist = {};\)\ \ \), "\n", 
    \(\(repeat = 200;\)\), "\n", 
    \(Print["\<actual = \>", actual]\), "\n", 
    \(\(count = 0;\)\), "\n", 
    \(Do[{count = 0, 
        Do[{Clear[x, y, z], x = Random[Real, {xmin, xmax}], 
            y = Random[Real, {ymin, ymax}], 
            z = Random[Real, {zmin, zmax}], \n\t\t\ \ \ \ \ \ If[z < f[x, y], 
              count = count + 
                  1]}, {numberofpoints}], \n\t\t\ \ \ \ \ \ estimate = 
          N[volumeofbox\ \ count/numberofpoints], error = estimate - actual, 
        AppendTo[errorlist, error]}, \n\t{repeat}]\), "\n", 
    \(\(percenterrorlist = errorlist/actual;\)\), "\n", 
    \(\(ListPlot[percenterrorlist, \ PlotJoined \[Rule] True];\)\)}], "Input"],

Cell["\<\
Look at the pattern of your errors. These are referred to as random errors, \
and they should be centered around 0 and should show no particular pattern. \
We can order them and look at that picture as well.\
\>", "Text"],

Cell[BoxData[{
    \(sortlist = Sort[percenterrorlist]\), "\n", 
    \(\(ListPlot[sortlist, \ PlotJoined \[Rule] True];\)\)}], "Input"],

Cell["\<\
You should run this again, increasing the number of repeats and/or the number \
of points to see what happens to the errors. Which will decrease your error \
size?\
\>", "Text"],

Cell["\<\
To analyze the nature of these errors,we will order them and display them in \
a bar chart.Be sure to execute the command for reading in the packages ONLY \
ONCE.\
\>", "Text"],

Cell[BoxData[{
    \(<< Statistics`DescriptiveStatistics`\n
    << Statistics`DataManipulation`\), "\n", 
    \(<< Graphics`Graphics`\)}], "Input",
  PageWidth->Infinity],

Cell["\<\
As you prepare the bar chart, you will probably have to adjust your intervals \
for the min and max and width. The min and max should be rounded down and up, \
respectively, to give \"nice\" values. The width is important in deciding how \
many bins you will have. The BinCount command determines how many data points \
fall in each bin and identification of the midpoints of those bins assists in \
the plot of the bar chart.
The LocationReport gives various measures of central tendency and the \
DispersionReport gives various measures of variation of the data.\
\>", "Text"],

Cell[BoxData[{
    \(\(min = \(- .30\);\)\), "\n", 
    \(\(max =  .30;\)\), "\n", 
    \(\(width =  .05;\)\), "\n", 
    \(frequencies = BinCounts[percenterrorlist, {min, max, width}]\), "\n", 
    \(midpts = 
      Table[min + width/2 + i\ width, {i, 0, 
          Length[frequencies] - 1}]\), "\n", 
    \(bc = 
      BarChart[Transpose[{frequencies, midpts}], 
        PlotRange \[Rule] All]\), "\n", 
    \(LocationReport[percenterrorlist]\), "\n", 
    \(DispersionReport[percenterrorlist]\)}], "Input"],

Cell["\<\
You should run this again with the number of times you repeat the method \
increased say from 200 to 500 and see what happens to the histogram. Does it \
become more bell-shaped?\
\>", "Text"]
}, Closed]],

Cell[CellGroupData[{

Cell["You Try It: Part II", "Section"],

Cell[TextData[{
  "You can extend the problem you began in the previous You Try It section by \
analyzing your ",
  StyleBox["errorlist",
    FontWeight->"Bold"],
  " when you repeat the Monte Carlo method with 500 points being generated \
each of 100 times."
}], "Text"],

Cell[BoxData[{\(Clear[x, y, z, errorlist, \ percenterrorlist]\), 
    RowBox[{
      RowBox[{"numberofpoints", "=", 
        StyleBox["500",
          FontColor->RGBColor[1, 0, 0]]}], ";", "\n", \(errorlist = {}\), 
      ";", "  ", "\n", \(repeat = 100\), ";", "\n", 
      \(Print["\<actual = \>", actual]\)}], 
  \(count = 0; \n
    Do[{count = 0, 
        Do[{Clear[x, y, z], x = Random[Real, {xmin, xmax}], 
            y = Random[Real, {ymin, ymax}], z = Random[Real, {zmin, zmax}], \n
            \t\t\ \ \ \ \ \ If[z < g[x, y], count = count + 1]}, {
            numberofpoints}], \n\t\t\ \ \ \ \ \ 
        estimate = N[volumeofbox\ \ count/numberofpoints], 
        error = estimate - actual, AppendTo[errorlist, error]}, \n\t{repeat}]
      \), \(percenterrorlist = errorlist/actual; \n
    ListPlot[percenterrorlist, \ PlotJoined \[Rule] True]; \)}], "Input"],

Cell[TextData[{
  "Are your errors centered around 0? If not, there is a mistake here.\nSee \
what happens to the range of your errors when you increase the ",
  StyleBox["numberofpoints",
    FontWeight->"Bold"],
  " command above."
}], "Text"]
}, Closed]]
}, Open  ]]
},
FrontEndVersion->"4.1 for Microsoft Windows",
ScreenRectangle->{{0, 1024}, {0, 695}},
WindowSize->{412, 433},
WindowMargins->{{5, Automatic}, {Automatic, 5}},
PrintingCopies->1,
PrintingPageRange->{Automatic, Automatic},
StyleDefinitions -> "Default.nb"
]

(*******************************************************************
Cached data follows.  If you edit this Notebook file directly, not
using Mathematica, you must remove the line containing CacheID at
the top of  the file.  The cache data will then be recreated when
you save this file from within Mathematica.
*******************************************************************)

(*CellTagsOutline
CellTagsIndex->{}
*)

(*CellTagsIndex
CellTagsIndex->{}
*)

(*NotebookFileOutline
Notebook[{
Cell[1705, 50, 362, 6, 109, "Text"],

Cell[CellGroupData[{
Cell[2092, 60, 123, 3, 390, "Title"],
Cell[2218, 65, 118, 3, 37, "Text"],

Cell[CellGroupData[{
Cell[2361, 72, 31, 0, 59, "Section"],
Cell[2395, 74, 143, 3, 52, "Text"],
Cell[2541, 79, 456, 7, 90, "Text"],

Cell[CellGroupData[{
Cell[3022, 90, 74, 1, 47, "Subsection"],
Cell[3099, 93, 1209, 34, 242, "Text"]
}, Closed]]
}, Closed]],

Cell[CellGroupData[{
Cell[4357, 133, 49, 0, 39, "Section"],
Cell[4409, 135, 1799, 35, 470, "Text"],
Cell[6211, 172, 627, 16, 313, "Input"],
Cell[6841, 190, 523, 11, 128, "Text"],
Cell[7367, 203, 56, 1, 30, "Input"],
Cell[7426, 206, 1237, 24, 610, "Input"],
Cell[8666, 232, 215, 4, 71, "Text"],
Cell[8884, 238, 1162, 23, 590, "Input"],
Cell[10049, 263, 124, 3, 52, "Text"]
}, Closed]],

Cell[CellGroupData[{
Cell[10210, 271, 37, 0, 39, "Section"],
Cell[10250, 273, 1056, 27, 280, "Text"],
Cell[11309, 302, 1087, 26, 290, "Input"],
Cell[12399, 330, 1018, 18, 610, "Input"]
}, Closed]],

Cell[CellGroupData[{
Cell[13454, 353, 42, 0, 39, "Section"],
Cell[13499, 355, 247, 4, 90, "Text"],
Cell[13749, 361, 824, 16, 450, "Input"],
Cell[14576, 379, 231, 4, 71, "Text"],
Cell[14810, 385, 135, 2, 50, "Input"],
Cell[14948, 389, 187, 4, 71, "Text"],
Cell[15138, 395, 186, 4, 71, "Text"],
Cell[15327, 401, 170, 4, 70, "Input"],
Cell[15500, 407, 588, 9, 166, "Text"],
Cell[16091, 418, 509, 12, 270, "Input"],
Cell[16603, 432, 202, 4, 71, "Text"]
}, Closed]],

Cell[CellGroupData[{
Cell[16842, 441, 38, 0, 39, "Section"],
Cell[16883, 443, 271, 7, 71, "Text"],
Cell[17157, 452, 871, 16, 450, "Input"],
Cell[18031, 470, 245, 6, 71, "Text"]
}, Closed]]
}, Open  ]]
}
]
*)



(*******************************************************************
End of Mathematica Notebook file.
*******************************************************************)