Thread: monte carlo methods

  1. #1
    Registered User
    Join Date
    Sep 2001
    Posts
    85

    monte carlo methods

    i have to modify this program to calculate the PI. I understand basicaly that the monte carlo methods are used to calculate numbers in a statistical simulation, but i don't know how to apply the methods to the problem. here is the program:



    // MONCARLO.CPP
    //
    // This program calculates the perimeter of an ellipse.
    // One quarter of the perimeter of an ellipse with the
    // larger radius 1 is equal to the so-called
    // "elliptic integral" E -- the area under the curve
    // y = sqrt(1 - (k*sin(x))^2) on the interval 0 <= x <= pi/2.
    // The parameter k describes the shape of the ellipse and is
    // usually expressed in terms of the so-called "modular angle":
    // k = sin(alpha), where alpha is the angle between the minor
    // axis and the line that connects its end to a focus of
    // the ellipse. For a circle, alpha = 0 and k = 0.
    // The integral is calculated by using the Monte Carlo method.
    //
    // Author: B. Random
    // Date: 2/29/1999
    //

    #include <iostream.h>
    #include <iomanip.h>
    #include <stdlib.h> // Declares rand() and RAND_MAX
    #include <math.h> // Declares sin(x)

    const double PI = 3.14159265358979323846;

    // Modular angle (in degrees and radians):
    const double alpha = 45;
    const double alpha_radians = alpha * PI / 180.;

    // k defines elongation of the ellipse:
    const double k = sin(alpha_radians);

    //************************************************** **************

    double MonteCarlo (long nPoints)

    // This function calculates the elliptic integral E.
    // It generates "nPoints" (pseudo)random points in the
    // rectangle { 0 <= x <= pi/2; 0 <= y <= 1 }, finds
    // the fraction of all points that fall under the curve,
    // and calculates the corresponding fraction of the area of
    // the rectangle.

    {
    double x, y, f;
    long n, count = 0;

    for (n = nPoints; n > 0; n--) {

    // Generate a random point:
    // 0 <= x <= pi/2
    // 0 <= y <= 1

    x = double(rand()) * (.5 * PI / RAND_MAX);
    y = double(rand()) / RAND_MAX;

    // Calculate 1 - (k*sin(x))^2

    f = k * sin(x);
    f *= f; // Square f
    f = 1 - f;

    // Increment the counter if the point (x,y) is on or under
    // the curve:

    if (y*y <= f) count++;
    // (more efficient than: if (y <= sqrt(f)) ...)
    }

    // count E
    // ------- = -----------------------------------------
    // nPoints area of rect {0 <= x <= pi/2; 0 <= y <= 1}

    return (double(count) / double(nPoints)) * (.5 * PI);
    }

    //************************************************** **************

    int main()

    {
    long n;
    double ellipsePerim;

    cout << "Enter the number of random points in Monte Carlo ==> ";
    cin >> n;
    ellipsePerim = 4. * MonteCarlo(n);

    cout << "The perimeter of the ellipse with a modular angle of "
    << alpha << " degrees\n"
    << " and a major axis of 2 units is approximately = "
    << setprecision(3) << ellipsePerim << endl;

    return 0;
    }

  2. #2
    Blank
    Join Date
    Aug 2001
    Posts
    1,034
    I'm not sure if this will help, but one way you could do it is
    instead of ellipse do a circle. Calculate the area
    by repeated taking a random number in a 1x1 square inscribed
    with a circle, and seeing if it's in the circle. Then it should be
    fairly easy to find PI because A = PI*r*r

  3. #3
    Registered User
    Join Date
    Sep 2001
    Posts
    85
    thats what i was thinking, so that means i need to make 0<=x and y <=1, right?

    so:

    x = double(rand()) / RAND_MAX;
    y = double(rand()) / RAND_MAX;


    but i don't understand what this means:
    f = k * sin(x);
    f *= f; // Square f
    f = 1 - f;

  4. #4
    Blank
    Join Date
    Aug 2001
    Posts
    1,034
    Basically it calculates E and then compares it to
    y to see if it's in the graph's area,
    but a speed up is realized so it calculates E^2 and
    compares to y^2.

    What you want to do with the circle is to see if it's in the inscribed circle.

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Generalising methods to handle multiple weapon types
    By Swarvy in forum Game Programming
    Replies: 2
    Last Post: 05-22-2009, 02:52 AM
  2. Class methods in header or source file?
    By TriKri in forum C++ Programming
    Replies: 13
    Last Post: 09-17-2007, 05:23 AM
  3. Turn Based Stradegy System Methods
    By TylerMoyer in forum Game Programming
    Replies: 2
    Last Post: 07-30-2007, 10:45 PM
  4. Lesson #5 - Methods
    By oval in forum C# Programming
    Replies: 1
    Last Post: 05-04-2006, 03:09 PM
  5. Replies: 8
    Last Post: 07-27-2003, 01:52 PM