Using a Monte Carlo Method for Definite Integration

Steven Keltner and Scott Sloan
University of Missouri-Rolla

Mathematica in Education
Volume 3 Number 1
Winter 1994
(c) TELOS/Springer-Verlag

Introduction

During our second semester Calculus class, we learned about many types of integration.  After learning about integrating by using antiderivatives, we explored numerical integration with the Trapezoid Method and Simpson's Rule.  Our professor mentioned that it would be possible to evaluate definite integrals by generating random points and determining whether or not they fell under the area of the function.  He said that a Mathematica  program could be created that would perform the operations.  We were very interested in exploring the idea of using a Monte Carlo method to  evaluate definite integrals.  Since we wanted to test the accuracy and precision of the method, it was necessary to design a program that would perform the procedure several times with varying numbers of random points for each calculation.

Programming a Monte Carlo Method of Integration

The program evaluates definite integrals by creating a "box" with dimensions defined by the limits of integration and the maximum and minimum values of the function on the selected interval.  The computer randomly generates the coordinates of a point within the box.  The function is evaluated at the randomly selected x value.  If the random point falls within the area of the integral, it is counted by Mathematica.  After performing the procedure for the specified number of points (iterations), the computer calculates the final result by multiplying the total area of the box by the percentage of points that fell within the bounds of the integral.  Of course, many functions have regions that fall below the line y=0.  If the program encounters a point that is within the area of the integral but below the x-axis, it subtracts one from the total number of points counted.
     One would expect that the more points generated for each evaluation of the integral, the more accurate the answer will be.  Since random numbers are used, however, there will always be variation from one determination to the next.  Because of the lack of precision, the program is designed to perform many calculations with a varying number of iterations.  Doing so makes it possible to "zero in on" the actual answer.
     The first part of the program defines the parameters:

f[x_]= Sin[x];
begint= -1;
endint= 1;
minit= 100;
maxit= 300;
stepsize= 100;
trial= 3;

The function is input as f[x], while the limits of integration are specified by begint and endint.  Since the program graphs the approximated value of the definite integral versus the number of iterations (points) used, the minimum and maximum number of iterations must be input as minit and maxit .  The distance between successive iteration values is determined by the stepsize.  The number of calculations performed at each value is specified by the trial statement. 	Thus, given the parameters above, the computer would use 100 random points in each of three separate determinations, 200 points each in three additional calculations, and 300 points for three final computations.
     The following lines determine the maximum and minimum value of the function on the specified interval.  The program uses these results as boundaries between which random y values are selected.  The values are also necessary to determine the area of the "box."

miny = Min[Table[N[f[x]],{x,begint,endint,0.01}]];
maxy = Max[Table[N[f[x]],{x,begint,endint,0.01}]];

The next part of the program  generates random coordinates and determines how each point should be counted.  If the random y value r is less than or equal to the function evaluated at the random x value s, it is counted.

Do[Do[Counter = 0;
Do[r = Random[Real, {miny,maxy}]; 
s = Random[Real,{begint,endint}];
If[r <= f[s] && r >= 0, Counter=Counter+1, Counter=Counter];
If[r >= f[s] && r < 0, Counter=Counter-1, Counter=Counter];
     {i,1,k}];

The portion of the program that follows determines the final result for each trial.  The number of points counted is divided by the total number of points, k, and multiplied by the area of the box.  The process is repeated for values of k ranging from minit to maxit.

AppendTo[v,N[(Counter/k)*((endint-begint)*(maxy-miny))]];
AppendTo[z,k],{j,1,trials}];
,{k,minit,maxit,stepsize}];

After transposing the data, the information is plotted:

data=Transpose[{z,v}];
ListPlot[data,PlotRange->All];

In order to test the accuracy of the results, we ran the program for a function that could be easily integrated.  The actual value of the definite integral was plotted along with the Monte Carlo results by replacing the two lines above with the following commands:

data=Transpose[{z,v}];
points=ListPlot[data,DisplayFunction->Identity,PlotRange->All];
value=NIntegrate[f[x],{x,begint,endint}];
actualval=Plot[value,{x,minit - 25, maxit},DisplayFunction->
     Identity];
Show[points,actualval,DisplayFunction->$DisplayFunction]

The parameters used in Graph 1 are as follows:
     f[x] = x^2 - 1
     begint = 0
     endint = 2
     minit = 1000
     maxit = 9000
     stepsize = 2000
     trials = 9

Graph 1.   Monte Carlo Integration of  x2 - 1

The graph demonstrates that the accuracy and precision of the results change as the number of points used in each calculation increases.  Generally, the more points that are used, the more accurate and precise the output will be.
     One of the most useful applications of the program is the evaluation of definite integrals that would otherwise be very difficult to calculate.  Graph 2 is a plot of the function Sin[1/x] on the interval -1 ² x ² 3.  Although Mathematica cannot integrate the function on this interval, results may be obtained using the Monte Carlo method.  Because the program uses the Min and Max commands, however,  it does not work well for intervals having an endpoint that is undefined. When this problem occurs, we recommend choosing an endpoint that is slightly larger or smaller than the original value (for example, select 0.001 instead of 0 for the function Sin[1/x] ).  Graph 3 is a plot of the results obtained for the definite integral of Sin[1/x] with -1 ² x ² 3.

Plot[Sin[1/x],{x,-1,3}]

Graph 2.   Sin[1/x]  for  -1 ² x ² 3

Graph 3.   Monte Carlo Integration of Sin[1/x] With -1 ² x ² 3

The Program

The following is a listing of the entire program.  It may be useful to insert

Print[N[(Counter/k)*((endint-begint)*(maxy-miny))]]

inside the second Do loop so that Mathematica displays the numerical result for each iteration as it is computed.

f[x_]=  ;    (*function = ??*) 
begint =  ;  (*begin interval*)
endint =  ;   (*end interval*)
minit =  ; (*smallest # iterations*)
maxit =  ; (*largest # iterations*)
stepsize =  ; (*step size*)
trials =  ;   (*repeats for same # iterations*)
z = {}; (*empty list*)
v = {}; (*empty list*)

miny = Min[Table[N[f[x]],{x,begint,endint,0.01}]];
maxy = Max[Table[N[f[x]],{x,begint,endint,0.01}]];

Do[Do[Counter = 0;
Do[r = Random[Real, {miny,maxy}]; 
s = Random[Real,{begint,endint}];
If[r <= f[s] && r >= 0, Counter=Counter+1, Counter=Counter];
If[r >= f[s] && r < 0, Counter=Counter-1, Counter=Counter],{i,1,k}];

AppendTo[v,N[(Counter/k)*((endint-begint)*(maxy-miny))]];
AppendTo[z,k],{j,1,trials}];
,{k,minit,maxit,stepsize}];

data=Transpose[{z,v}];
ListPlot[data,PlotRange->All];

The Authors

Steven Keltner
Rt 3 Box 515
Willard, MO  65781
(417) 742-2311

Scott Sloan
Rt 2 Box 146
Plattsburg, MO  64477
(816) 539-2007