Computing integral using Monte Carlo method

[gfd43tg]

Homework Statement

Complete the following MATLAB function that computes the integral $\int_{-2}^{5} f(x)dx$ of the function
$f(x)$ plotted in the figure below, after N trials using a Monte Carlo method. Take all the required
numerical inputs from the figure shown below.
Hint: Generate a set of $N$ uniformly distributed random points $(x, y)$ within a rectangular area
and determine the fraction of them that fall within the area of interest (i.e., under $f(x)$).

function val = mcIntegral(fh,N)
% Inputs: fh: the function handle to vectorized function f(x)
% -2, 5: the lower and upper integration bounds
% N: the number of trials
% Output: val: the calculated value of the integral

A =_______________________________________________________;% enclosing area


x = rand( )
___________________________________________________________________


y = rand( )
___________________________________________________________________


success =_________________________________________________________;
val = A * success / N ;

Homework Equations

The Attempt at a Solution

I have never heard of Monte Carlo, but this was given on a previous exam. I am wondering if one could reasonably answer this question without ever hearing about the method, and how I will be able to do it?? It obviously has something to do with probability. My guess is A is supposed to be a rectangle with a width of $5- (-2) = 7$ by height of the max value of the function, but I don't know how to get that maximum. I also don't know what to do with success, or if x and y are okay.

function val = mcIntegral(fh,N)
A = ? [7 max fh(x)]
x = rand(N,1);
y = rand(N,1);
success = (x <= 5) && (x >= -2) &&  (y <= fh(x))
val = A*success/N

 

[/AH]

Monte Carlo is about the probability that a random point (x,y) will lie within your desired area. So the desired area is a fraction of the total area in which you evaluate points. Do you have an idea what those two area's might be? And if so, how do you determine if a point is within your desired area?

 

[gfd43tg]

I figure the total area is a rectangle that has a width of 7 and a value of the function at a specific x. I determine that it is within the desired area with my conditional statements

 

[AlephZero]

You don't need to know the max and min values of f(x) beforehand, if you do the steps in the right order.

1. Generate a random set of x values.
2. Find the corresponding values of f(x).
3. Set fmax and fmin to the maximum and minimum of the values from step 2.
4. Generate a random set of y values between fmax and fmin.
5. Count the number of times the point $(x_i,y_i)$ is below the point $(x_i,f(x_i)$.
6. Find the approximate value of the integral.

Note step 6 isn't quite straightforward. If fmin > 0, the integral also needs to include the area of the rectangle between y = 0 and y = fmin. You also need to deal with the situations when fmin < 0 and fmax > 0, and when fmax < 0.

 

[gfd43tg]

Yes, but I have to follow the format given in the problem statement, so how do I do that?

 

[pbuk]

Look at the plot of the function in the question, it shows clearly the intended rectangle.

x = rand(N,1) generates a vector of random numbers in the range (0,1) but you want the range (-2,5).
y = rand(N,1) generates a vector of random numbers in the range (0,1) but you want the range [look at the plot].
You are close with success, but you are currently setting it to an array of logical 0/1 values: what should it be instead?

 

[gfd43tg]

MrAnchovy,

I am inexperienced with Monte Carlo methods. This ''obvious'' rectangle is a complete mystery to me, unless you mean the entire graph is the intended rectangle. As for y, what is the range of the plot? Just the maximum value of the function? Does it matter if my rectangle is out of bounds of the plot?

so x = -2 + 7*rand(N,1) ?? But for y I still don't know

 

[/AH]

The integral of the function f(x) on range(-2 to 5) is the area enclosed by the the function and the x-axis, as you most likely are familiar with. A hint to the answer is in the final rule in the code:
val = A * success / N;
Your answer is a total area times the fraction of "succesfull" points.

So on one hand you have a "total area" on the other hand you have your answer which is the area under the function f(x). So as MrAnchovy point out as well, the most obvious choice for your total area would simply be the entire rectangular enclosure of the graph.

The success rate you can determine as Alephzero said, by evaluating your random point by:
$(x_i,y_i)<(x_i,f(x_i))$

EDIT: Note, your choice for the total area arbitrary as long as it fully encloses your desired answer. Selection of the area is based on speed. A larger area (i.e. y-range(0,100) means more points to evaluate). You can in theory also change your x-range, but that would require evaluating in x-direction as well which is not advantageous.

 

[gfd43tg]

How is this? I know it's not right, I think success is not quite there yet

function val = mcIntegral(fh,N)
A = 7*10;
x = -2 + 7*rand(N,1);
y = 10*rand(N,1);
success = numel((x <= 5) && (x >= -2) &&  (y <= fh(x)));
val = A*success/N;

 

[pbuk]

Close. numel will return N so you want a different function to sum the logical 1 values. It should then work, although your answer should be penalised for redundant code; remove two inequalities.

 

[gfd43tg]

function val = mcIntegral(fh,N)
A = 7*10;
x = -2 + 7*rand(N,1);
y = 10*rand(N,1);
success = sum(y <= fh(x));
val = A*success/N;

That wasn't all that bad after all! But if I had never seen Monte Carlo and it was on the exam I would have been a goner I think.

 

[pbuk]

The term "Monte Carlo" is applied to many methods for numerical solution of various types of problem where there is an element of random selection rather than calculation of sample points, step size, covering grids etc. Monte Carlo methods are particularly suited to problems involving functions with an unknown periodic element where the randomness tends to avoid problems with bias in finite element selection.

What course are you doing? I'm surprised you are preparing from what appear to be computational maths exam papers without having come across the term.

 

[gfd43tg]

Introduction to programming in Matlab