assassin_696 Posted February 17, 2011 Share Posted February 17, 2011 I'm trying to write a program that calculates the value of ln(2) using monte carlo methods in C++. Obviously ln(2) can be calculated from the integral of 1/x between the limits 2 and 1, but I'm having some trouble implementing that into code. So far I've got: //Program uses Monte Carlo (random number generator) methods //to calculate a value for ln(2) #include <iostream> #include <cstdlib> #include <cstdio> using namespace std; #define ranf() ?????????? int main(int argc, char* argv[]) { int outcome, N=10, count_in=0, seed=30; double fraction_in; if(argc>1) { sscanf( argv[1], "%d", &N ) ; // puts first command line argument in N } if(argc>2) { sscanf( argv[2], "%d", &seed ); // puts the 2nd argument in seed } cout << " # " << argv[0] << " N= " << N << " seed= " << seed << endl; srand(seed); // initialises the random number generator for(int n=1; n<=N; n++) { double x =ranf(); // ranf() needs to return a real number in [0,1] double y =ranf(); outcome = (x*y < 1.0 ) ? 1 :0; count_in += outcome; fraction_in = static_cast<double>(count_in)/n; //casts integer variables for correct division cout << "Location" << outcome << "\t" << x << "\t" << y << "\t" << count_in << "\t" << n << "\t" << fraction_in <<endl; } return 0; } The trouble I'm having is putting the limits of the integration in, which I think is done in the ranf() function where I put a load of question marks. Looking at a similar example for calculating the value of pi by testing whether or not a point is in the unit circle, ranf is #define ranf() \ ((double)rand()/(1.0+(double)0x7fffffff)) Which I think is where the limits information is obtained, but I'm not sure. I also want to make sure that my outcome expression, i.e. outcome = (x*y < 1.0 ) ? 1 :0; Is correct, since that's where the integral equation information is contained, obviously if I screwed that up it'd be calculating the wrong integral. Any help would be appreciated, I'm pretty new to all this so not sure how easy that code is to read for more experienced people but obviously I can clarify anything. "Da mihi castitatem et continentam, sed noli modo" Link to comment Share on other sites More sharing options...
Markup Posted February 18, 2011 Share Posted February 18, 2011 (edited) Updating, ignore http://www.cplusplus.com/reference/std/complex/log/http://www.cplusplus.com/reference/clibrary/cstdlib/rand/ The second link being the most relevant. It shows you how to setup a seeded random generator. It also explains how to limit the range using modulo. replaced! Enjoy Notes: Observe indentation, although opinion based ;pThe declaration of x and y are outside of the for loop to stop performance impact of declaring the variable on each iteration of the loop. outcome = (x*y < 1.0 ) ? 1 : 0; If condition is true, return the first value(left side of : ), if false return the second value(right side of : ) Edited February 19, 2011 by Markup Link to comment Share on other sites More sharing options...
Veiva Posted February 18, 2011 Share Posted February 18, 2011 In order to generate a number in the range of [0, 1], the usual best method is to divide the return value of rand by the value of RAND_MAX. Like so: #define randf() (rand() / (double)RAND_MAX) That will assuredly give you a value between 0 and 1. You can verify the math yourself using a set using an example implementation of rand that has a maximum return value of 255. For any value between 0 and 255, divide it by 255. You'll get 0, 0.004, 0.008, ..., 1. Link to comment Share on other sites More sharing options...
dsavi Posted February 19, 2011 Share Posted February 19, 2011 Wow I wish I understood this. Link to comment Share on other sites More sharing options...
Markup Posted February 19, 2011 Share Posted February 19, 2011 Wow I wish I understood this. I understand the monte carlo concept, just unsure on how to implement the finer details. In layman's terms,You place a box around the area of a graph you want to integrate, then using a random sampling method, you test points in the box. You keep track of how many points tested for and how many are under the line. Using this data you can approximate areas or expressions. eg/An illustration of Monte Carlo integration. In this example, the domain D is the inner circle and the domain E is the square. Because the square's area can be easily calculated, the area of the circle can be estimated by the ratio (0.8) of the points inside the circle (40) to the total number of points (50), yielding an approximation for pi/4 approx 0.8 (PI*(d^2))/4 = (40/50)*(d^2) => PI/4 ~ (40/50) The calculated approx area being, (40/50)*(d^2) = 3.2 Link to comment Share on other sites More sharing options...
Markup Posted February 19, 2011 Share Posted February 19, 2011 Think I got it //Program uses Monte Carlo (random number generator) methods //to calculate a value for ln(2) #include <iostream> #include <stdio.h> #include <stdlib.h> #include <time.h> #define myRand() ((double)rand() / (double)RAND_MAX); int main(int argc, char* argv[]) { int seed = (int)time(NULL); int n=100000; double x; double y; int r = 0; double ratio = 0.0; srand(seed); for(int i = 0; i < n; i++) { x = myRand(); y = myRand(); x += 1.0; r += (y*x <= 1.0) ? 1 : 0; //uncomment below to enter the matrix! //ratio = (double)r/(double)i; //std::cout << "Iteration: " << i << "\tArea: " << ratio*2.0 << "\tRatio: " << ratio << "\tx: " << x << "\ty: " << y << std::endl; } ratio = (double)r/(double)i; std::cout << "Ratio: " << ratio << "\tArea: "<< ratio*2.0 << std::endl; return 0; } Link to comment Share on other sites More sharing options...
Recommended Posts
Create an account or sign in to comment
You need to be a member in order to leave a comment
Create an account
Sign up for a new account in our community. It's easy!
Register a new accountSign in
Already have an account? Sign in here.
Sign In Now