Take Your Chances: Try the Monte Carlo Technique for Numerical Integration in Three Dimensions

Introduction

OBJECTIVE: Learn the Monte Carlo method for approximating a multiple integral that cannot be integrated symbolically.

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.

Technology Guidelines

NOTE:  If you have just finished a module, restart Mathematica or Quit Kernel under the Evaluation pull-down menu before executing a new module.
TO OPEN CELLS, click on the arrow to the left of the cell or put your cursor on the right cell bracket and double click. To open all cells, highlight the bracket on the far right and select Grouping under the Cell pull-down menu and select Open All Subgroups.
TO STOP AN EXECUTION
Select the Evaluation pull-down menu and click on Abort Evaluation.
ORDER OF EXECUTION
Execute cells in the order given. Do not skip any Input cells within a given notebook.
SAVING NOTEBOOKS
You can save anytime to any directory you choose, and it is wise to save often. However, before you do your final save to reduce the size of your file, it is a good idea to delete all your output by selecting the Delete All Output selection under the Cell pull-down menu.
EXPERIENCING MAJOR PROBLEMS
Save if appropriate, then shut down Mathematica and start it up again.

Part I: The Monte Carlo Method

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 Mathematica'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.

In utilizing the Random command in Mathematica, the default Random[] 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 x-y 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

Clear[x,y,z,f]
f[x_,y_]:=(E^(-2 (x^2)-(y^2)))
xminf=0;
xmaxf=1.5`;
yminf=0;
ymaxf=2;
zminf=0;
zmaxf=1;
pf=Plot3D[f[x,y],{x,xminf,xmaxf},{y,yminf,ymaxf},PlotRange->All]
actualf=NIntegrate[f[x,y],{x,xminf,xmaxf},{y,yminf,ymaxf}];
Print["actual volume = ",actualf]

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 listin are under the surface in the volume we are computing and they will be plotted in red. The points in listout are not under the surface and they will be plotted in black.

numberofpoints=1000;
volumeofboxf=(xmaxf-xminf) (ymaxf-yminf) (zmaxf-zminf);
count=0;
listout={};
listin={};
Do[{Clear[x,y,z],SeedRandom[],x=RandomReal[{xminf,xmaxf}],y=RandomReal[{yminf,ymaxf}],z=RandomReal[{zminf,zmaxf}],If[z≤f[x,y],{count=count+1,AppendTo[listin,{x,y,z}]},AppendTo[listout,{x,y,z}]]},{numberofpoints}]
pout=ListPointPlot3D[listout,DisplayFunction->Identity];
pin=ListPointPlot3D[listin,PlotStyle->RGBColor[1,0,0],DisplayFunction->Identity];
Show[pin,pout,pf,DisplayFunction->$DisplayFunction]
Print["actual = ",actualf]
estimate=N[((volumeofboxf count)/numberofpoints)];
Print["estimate = ",estimate]
error=estimate-actualf;
Print["error = ",error]
relativeerror=(error/actualf);
Print["relative error = ",relativeerror]

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.

Clear[x,y,z]
numberofpoints=3000;
count=0;
listout={};
listin={};
Do[{Clear[x,y,z],x=RandomReal[{xminf,xmaxf}],y=RandomReal[{yminf,ymaxf}],z=RandomReal[{zminf,zmaxf}],If[z≤f[x,y],{count=count+1,AppendTo[listin,{x,y,z}]},AppendTo[listout,{x,y,z}]]},{numberofpoints}]
pout=ListPointPlot3D[listout,DisplayFunction->Identity];
pin=ListPointPlot3D[listin,PlotStyle->RGBColor[1,0,0],DisplayFunction->Identity];
Show[pin,pout,pf,DisplayFunction->$DisplayFunction]
Print["actual = ",actualf]
estimate=N[((volumeofboxf count)/numberofpoints)];
Print["estimate = ",estimate]
error=estimate-actualf;
Print["error = ",error]
relativeerror=(error/actualf);
Print["relative error = ",relativeerror]

Does this method seem very accurate?  We will analyze the error involved in such an approximation.

You Try It: Part I

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 are your function and minimum and maximum values for each of your variables. These are all in red (leave your zmin at 0) . Then, simply execute the rest of the commands. Note that we are calling the function in the You Try It sections g(x,y) instead of f(x,y).

You 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?

We will begin by plotting the function and letting Mathematica estimate the volume between the surface and the xy-plane over the domain specified. Be sure to select a function that stays above the xy-plane.

Clear[x,y,z,g]
g[x_,y_]=((Sin[x-y])^2);
xming=-2;
xmaxg=2;
yming=0;
ymaxg=3;
zming=0;
zmaxg=1;
pfg=Plot3D[g[x,y],{x,xming,xmaxg},{y,yming,ymaxg},PlotRange->All,PlotPoints->40]
actualg=NIntegrate[g[x,y],{x,xming,xmaxg},{y,yming,ymaxg}];
Print["actual volume = ",actualg]

numberofpoints=2000;
volumeofboxg=(xmaxg-xming) (ymaxg-yming) (zmaxg-zming);
count=0;
listout={};
listin={};
Do[{Clear[x,y,z],x=RandomReal[{xming,xmaxg}],y=RandomReal[{yming,ymaxg}],z=RandomReal[{zming,zmaxg}],If[z≤g[x,y],{count=count+1,AppendTo[listin,{x,y,z}]},AppendTo[listout,{x,y,z}]]},{numberofpoints}]
pout=ListPointPlot3D[listout,DisplayFunction->Identity];
pin=ListPointPlot3D[listin,PlotStyle->RGBColor[1,0,0],DisplayFunction->Identity];
Show[pin,pout,pfg,DisplayFunction->$DisplayFunction]
Print["actual = ",actualg]
estimate=N[((volumeofboxg count)/numberofpoints)];
Print["estimate = ",estimate]
error=estimate-actualg;
Print["error = ",error]
relativeerror=(error/actualg);
Print["relative error = ",relativeerror]

Part II: Error Analysis

In this part, we follow the procedure above with 1000 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.

Clear[x,y,z]
numberofpoints=1000;
errorlist={};
repeat=200;
Print["actual = ",actualf]
count=0;
Do[{count=0,Do[{Clear[x,y,z],x=RandomReal[{xminf,xmaxf}],y=RandomReal[{yminf,ymaxf}],z=RandomReal[{zminf,zmaxf}],If[z≤f[x,y],count=count+1]},{numberofpoints}],estimate=N[((volumeofboxf count)/numberofpoints)],error=estimate-actualf,AppendTo[errorlist,error]},{repeat}]
Clear[relativeerrorlist]
relativeerrorlist=(errorlist/actualf);
ListPlot[relativeerrorlist,Joined->True]

Look at the pattern of your errors. These are referred to as random errors, and they will usually be centered around 0 and should show no particular pattern. We can order them and look at that picture as well.

sortlist=Sort[relativeerrorlist]
ListPlot[sortlist,Joined->True]

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?

You Try It: Part II

You can extend the problem you began in the previous You Try It section by analyzing your errorlist when you repeat the Monte Carlo method with 1000 points being generated each of 200 times.

Clear[x,y,z,errorlist,percenterrorlist]
numberofpoints=1000;
errorlist={};
repeat=200;
Print["actual = ",actualg]
count=0;
Do[{count=0,Do[{Clear[x,y,z],x=RandomReal[{xming,xmaxg}],y=RandomReal[{yming,ymaxg}],z=RandomReal[{zming,zmaxg}],If[z≤g[x,y],count=count+1]},{numberofpoints}],estimate=N[((volumeofboxg count)/numberofpoints)],error=estimate-actualg,AppendTo[errorlist,error]},{repeat}]
Clear[relativeerrorlist]
relativeerrorlist=(errorlist/actualg);ListPlot[relativeerrorlist,Joined->True]

Are your errors centered around 0? If not, there might be a mistake here.
See what happens to the range of your errors when you increase the numberofpoints command above.