1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Computing integral using Monte Carlo method

  1. Aug 13, 2014 #1

    Maylis

    User Avatar
    Gold Member

    1. The problem statement, all variables and given/known data
    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)##).

    Code (Text):
    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 ;
    2. Relevant equations



    3. 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.

    Code (Text):
    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
     

    Attached Files:

    Last edited: Aug 13, 2014
  2. jcsd
  3. Aug 14, 2014 #2

    /AH

    User Avatar

    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?
     
  4. Aug 14, 2014 #3

    Maylis

    User Avatar
    Gold Member

    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
     
  5. Aug 14, 2014 #4

    AlephZero

    User Avatar
    Science Advisor
    Homework Helper

    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.
     
  6. Aug 14, 2014 #5

    Maylis

    User Avatar
    Gold Member

    Yes, but I have to follow the format given in the problem statement, so how do I do that?
     
  7. Aug 14, 2014 #6
    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?
     
  8. Aug 14, 2014 #7

    Maylis

    User Avatar
    Gold Member

    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
     
    Last edited: Aug 14, 2014
  9. Aug 15, 2014 #8

    /AH

    User Avatar

    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:
    [itex](x_i,y_i)<(x_i,f(x_i))[/itex]

    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.
     
    Last edited: Aug 15, 2014
  10. Aug 15, 2014 #9

    Maylis

    User Avatar
    Gold Member

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

    Code (Text):
    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;
     
  11. Aug 15, 2014 #10
    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.
     
  12. Aug 15, 2014 #11

    Maylis

    User Avatar
    Gold Member

    Code (Text):
    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.
     
    Last edited: Aug 15, 2014
  13. Aug 15, 2014 #12
    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.
     
  14. Aug 15, 2014 #13

    Maylis

    User Avatar
    Gold Member

    Introduction to programming in Matlab
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Computing integral using Monte Carlo method
  1. Monte-Carlo Method (Replies: 2)

Loading...