topic:  The Monte Carlo method and multiple integration
description:  discusses the Monte Carlo method as a numerical method for approximating the values of multiple integrals.  Contains an application to weapon effectiveness.
file name:  montecarlo.nb
directions:  Read the text cells for background information.  Execute any existing input cells as you come to them.  Answer questions or work exercises as you come to them.  Modify the existing input cells as needed to work exercises requiring the use of Mathematica.

Computer Assignment #4
The Monte Carlo Method

Special thanks to Calculus&Mathematica (Addison-Wesley, 1994) for their support of these assignments.  Each is very much in the style of C&M, and some problems and commands are adaptations, with permission, of C&M problems and commands.

Why should I have to do this computer assignment anyway?

Are you still asking this?!
	The Monte Carlo method is an actual method used by scientists and engineers to gather statistical data on reliability and effectiveness of certain devices or processes.  Monte Carlo simulation is particularly useful for real-world systems involving random or near-random behavior.  Its great advantage is that Monte Carlo simulation can be performed on a computer in a relatively short amount of time, whereas actual field tests might require vast amounts of time, personnel, and equipment.

Introduction

In your last course, we talked about how to overapproximate and underapproximate the area between two curves by using Riemann sums.  In this computer assignment, we will talk about how to approximate both the area between two curves and the volume between two surfaces by using the Monte Carlo method.  Then, we will investigate an application of the Monte Carlo method.

Part 1 -- The Area of a Region

Suppose we wanted to APPROXIMATE the area of the region bounded above by
y = 2 - x^2 and below by y = x^4.  (Of course, we know how to get the exact answer...but let's pretend we don't.)   Execute the following cell to define the functions and to see the region:

Clear[f,g,x];
f[x_] = 2 - x^2;
g[x_] = x^4;
regionplot = Plot[{f[x],g[x]},{x,-1,1},
			 AspectRatio->Automatic,
			 AxesLabel->{"x","z"}];

We COULD approximate the shape of the region with a bunch of little rectangles, and then sum up the areas of the rectangles.  This is the basic idea behind a Riemann sum!  Alternatively, we could execute the following commands (ignore any "possible spelling error" messages):

NumberOfDarts = 500;
NumberOfHits = 0;
PointList = {};

xmin = -1;
xmax = 1;
zmin = 0;
zmax = 2;
RecArea = (xmax - xmin)*(zmax - zmin);

Randomize[];
Do[x = Random[Real,{xmin,xmax}];
   z = Random[Real,{zmin,zmax}];
   AppendTo[PointList,{x,z}];
   If[z >= g[x] && z <= f[x],
		NumberOfHits = NumberOfHits + 1],
   {i,1,NumberOfDarts}];

RegionArea = N[NumberOfHits/NumberOfDarts * RecArea];

Print["There were ",NumberOfHits," hits out of ",
	NumberOfDarts," test points."];
Print["The target rectangle has area ",RecArea,"."];
Print["The approximate area of the region is ",
	RegionArea,"."];
Clear[x,z];

Here is what the above algorithm does:
	The algorithm sets the variable NumberOfDarts, representing the total number of random points, at 500 and sets the variable NumberOfHits, representing the number of points which land within the region, initially at 0.  The algorithm also sets the dimensions of the target rectangle and calculates its area.
	The "Do" loop tells Mathematica to select NumberOfDarts random points {x,z}within the target rectangle and to check whether or not each is also within the region of interest.  For each one that is, Mathematica increments the counter NumberOfHits by 1.
	Finally, Mathematica approximates the area of the region and stores it in RegionArea.
	The following commands give a visual look at your Monte Carlo simulation above.  (The variable PointList and the AppendTo command are not really part of the Monte Carlo algorithm...those are included so that we can display the random points in the next command.)

darts = ListPlot[PointList,
		AspectRatio->Automatic,
		PlotStyle->PointSize[.015],
		DisplayFunction->Identity];
Show[regionplot,darts,DisplayFunction->$DisplayFunction];

a.  Suppose that you executed the Monte Carlo algorithm above and obtained the result that the approximate area of the region was 3.184.  Exactly how many of the 500 points would have landed within the region?

<<ANSWER QUESTION IN YOUR COPY OF THE FILE COMP4ANS.MA>>

b.  Execute the above algorithm 5 times.  Record the results of your trials below:

<<ANSWER QUESTION IN YOUR COPY OF THE FILE COMP4ANS.MA>>

c.  Based on these five trials, give a good estimate for the area of the region.  EXPLAIN how you arrive at this estimate.

<<ANSWER QUESTION IN YOUR COPY OF THE FILE COMP4ANS.MA>>

d.  Use the Integrate[] command to find the exact area of the region.

<<ANSWER QUESTION IN YOUR COPY OF THE FILE COMP4ANS.MA>>

e.  Now, set zmax = 30 in the above algorithm (instead of its current value of 2).  Run the algorithm, and report the estimated area of the region.

<<ANSWER QUESTION IN YOUR COPY OF THE FILE COMP4ANS.MA>>

f.  Is the estimate of Part e a better or a worse approximation than the estimates of
Part b?  What do we need to consider in setting up the Monte Carlo algorithm to find the area of a region?  EXPLAIN your reasoning.

<<ANSWER QUESTION IN YOUR COPY OF THE FILE COMP4ANS.MA>>

Part 2 -- The Volume of a Solid

The following plot depicts an inverted paraboloid and a saddle surface, plotted over the square -2 < x < 2, -2 < y < 2:

Clear[f,g,x,y];
f[x_,y_] = 15 - x^2 - y^2;
g[x_,y_] = x y;
fplot = Plot3D[f[x,y],{x,-2,2},{y,-2,2},
			DisplayFunction->Identity];
gplot = Plot3D[g[x,y],{x,-2,2},{y,-2,2},
			DisplayFunction->Identity];
Show[fplot,gplot,
	ViewPoint->{1.605,-2.964,0.297},
	DisplayFunction->$DisplayFunction];

Let's use the Monte Carlo method to approximate the volume between the two surfaces over this square.

a.  MODIFY the algorithm of Part 1 above to approximate the volume of the solid.  Realize that you now want to pick random test points in 3-space.

<<ANSWER QUESTION IN YOUR COPY OF THE FILE COMP4ANS.MA>>

b.  Use Mathematica to set up and evaluate an iterated double integral which gives the exact volume.

<<ANSWER QUESTION IN YOUR COPY OF THE FILE COMP4ANS.MA>>

Part 3 -- Why the Monte Carlo Method?

(Thanks to Capt James Wisnowski and Maj James Simpson for suggesting this problem and for providing technical assistance!)

Monte Carlo simulation is often used in the statistical analysis of weapon effectiveness.  This problem gives a bit of the flavor without getting too involved in gory details.
	Suppose we are testing the effectiveness of the Sensor Fuzed Weapon (SFW), used to neutralize enemy tanks.  The SFW consists of a main housing, called the Tactical Munitions Dispenser (TMD), which is released from an F-16.  The TMD deploys 10 submunitions, called Skeet Delivery Vehicles (SDVs), which each descend over a target area by parachute.  At a certain altitude, each SDV releases four explosive devices, called skeets.  The initial position of a skeet is essentially random.  As a skeet descends, it looks for the heat signature of a tank below it.  Once a heat signature is found, the skeet detonates.
	For simplicity in this problem, we will assume the skeets detonate at random positions.  Thus, our probabilities will be worse than the true probabilities.
	Suppose the target area is defined by the rectangle below, which contains three simulated enemy tanks:

drawtank[cx_,cy_] :=
	Graphics[{
	Line[{{cx+.16,cy+.01},
		{cx+.16,cy+.08},{cx-.16,cy+.08},
		{cx-.16,cy-.08},{cx+.16,cy-.08},
		{cx+.16,cy-.01}}],
	Circle[{cx,cy},.05],
	Line[{{cx+.05,cy-.01},{cx+.20,cy-.01},
		{cx+.20,cy+.01},{cx+.05,cy+.01}}]}];
tanks = Show[
	Graphics[Line[{{0,0},{2,0},{2,1},{0,1},{0,0}}]],
	drawtank[.4,.7],
	drawtank[.9,.3],
	drawtank[1.4,.75],
	AspectRatio->Automatic];

Now, we run the simulation.  The Monte Carlo simulation scatters 40 random points, each representing a skeet detonation, across the target area.  After the simulation, we count up the number of tanks hit by points.  If a particular tank contains a point, we consider the corresponding tank to be destroyed.

NumberOfSkeets = 40;
PointList = {};
xmin = 0;
xmax = 2;
ymin = 0;
ymax = 1;
Do[x = Random[Real,{xmin,xmax}];
   y = Random[Real,{ymin,ymax}];
   AppendTo[PointList,{x,y}],
   {i,1,NumberOfSkeets}];
hits = ListPlot[PointList,
		AspectRatio->Automatic,
		PlotStyle->PointSize[.015],
		DisplayFunction->Identity];
Show[tanks,hits,DisplayFunction->$DisplayFunction];
Clear[x,y];

a.  Run the above Monte Carlo simulation 10 times.  Count and record the number of tanks destroyed for each trial.

<<ANSWER QUESTION IN YOUR COPY OF THE FILE COMP4ANS.MA>>

b.  Based on your 10 trials, what is the probability that a single SFW with its 40 skeets will take out all 3 tanks on a given pass?

<<ANSWER QUESTION IN YOUR COPY OF THE FILE COMP4ANS.MA>>

c.  How would you use the Monte Carlo method to estimate the probability that a SINGLE skeet detonating at a random point within the rectangular target area hits one of the three tanks?

<<ANSWER QUESTION IN YOUR COPY OF THE FILE COMP4ANS.MA>>

During an actual deployment of the SFW, it is common for some of the skeets to fail, or for entire SDVs to fail.  A true Monte Carlo simulation factors in the probability of a component failing, so that generally there are not a full 40 detonations displayed.