## Calculating pi

## Calculating pi

(OP)

Hello everyone, I'm trying to calculate pi by following an idea on Wikipedia

http://en.wikipedia.org/wiki/Monte_Carlo_method

In practice I generate some random points (x,y) in a square (side=1) and counting how many points are in the square and how many are in the circle (radius=1) remembering that I'm considering a quarter of circle, not a full circle.

The problem is that I set an error threshold of 10^(-8) but the program stops after 120217 steps, with an error of 7*10^(-6).

How can I improve my program?

Thank you!

P.S. sorry for my bad English, I'm not mother tongue.

http://en.wikipedia.org/wiki/Monte_Carlo_method

In practice I generate some random points (x,y) in a square (side=1) and counting how many points are in the square and how many are in the circle (radius=1) remembering that I'm considering a quarter of circle, not a full circle.

The problem is that I set an error threshold of 10^(-8) but the program stops after 120217 steps, with an error of 7*10^(-6).

How can I improve my program?

Thank you!

P.S. sorry for my bad English, I'm not mother tongue.

#### CODE --> Fortran

program pi_montecarlo implicit none double precision::x,y,error,pi_calculated integer::i,n,square,circle double precision, parameter::pi=3.141592653589793238462643383279 square=0 circle=0 error=1. n=0 do while (abs(error)>1e-10) x=ran() y=ran() if (sqrt(x**2.+y**2.)<=1.) then circle=circle+1 square=square+1 else square=square+1 end if pi_calculated=4.*(float(circle)/float(square)) error=(pi_calculated-pi) write(*,*)"steps=",n,"pi=",pi_calculated,"error=",error n=n+1 end do end program

## RE: Calculating pi

## RE: Calculating pi

http://de.wikipedia.org/wiki/Kreiszahl#Statistisch...

and tried this

## CODE

Output:

## CODE

## RE: Calculating pi

double precision, parameter::pi=3.141592653589793238462643383279

does not provide the expected result : print pi to see the obtained value. You should be surprised.

The right declaration is :

double precision, parameter::pi=3.141592653589793238d0

Taking into account that

pi_calculated=4.*(float(circle)/float(square))

even if pi_calculated is declared "double precision" the right member is computed in single precision !

The right formula is :

pi_calculated=4.d0*(dble(circle)/dble(square))

At last, instead of rand(), you should use the standard subroutine random_number :

call random_number(x)

call random_number(y)

With that routine, you are sure that the return value is double precision (as the type of the argument). With rand() which is not standard, I don't know !

François Jacq

## RE: Calculating pi

Both of your answers were very helpful!

Thanks for the code, it's really useful to see another way to answer the problem.

And I didn't know I wasn't using double precision numbers for all my values. You've just opened a new world for me!

Thank you.

## RE: Calculating pi

n = 2147483647is the maximum value of the default integer type.If we need bigger n to achieve better approximation, we need to use other kind of integer.

## RE: Calculating pi

## RE: Calculating pi

http://gcc.gnu.org/onlinedocs/gfortran/SELECTED_00...

## RE: Calculating pi

## RE: Calculating pi

## CODE

## CODE

## CODE

It seems that the probability works slow...

Interesting would be the question how large should be the number of drops to get the approximation of PI with given accuracy, e.g. of 1e-10 (with probability of 99%)

## RE: Calculating pi

Although I have heard about Monte Carlo method, but until now I have never used it.

## RE: Calculating pi

=> to divide by 2 the distance to the right value, you need to multiply by 4 the number of attempts.

François Jacq

## RE: Calculating pi

I found this picture on Wikipedia that shows the percent error as a function of the number of drops, and it says that overall the program generated 1,37 billion drops, which is similar to the results you got with n=2147483647, as 100*(5,7E-5/pi)=0,0018%

http://upload.wikimedia.org/wikipedia/commons/7/7d...

Ohh and thanks to both of you for answering to my question; I'm glad it interested you too.

I tried my old program still using ran() and changing float in dble

## CODE --> Fortran

and got this

## CODE -->

but maybe it's just a case that it reached that precision in just a few steps, like the zero points in the graph on Wikipedia I posted, that can occur also for "very few" drops.

## RE: Calculating pi

## CODE

With n = 4678343 my code only computes this poor approximation

## CODE

## RE: Calculating pi

## CODE

## CODE

But when I call the function more times with several n - e.g. like this:

## CODE

## CODE

`drops_inside = 792 drops_total = 1000 n = 1000 pi = 3.1680000000000001 Err = 2.6407346410207033E-002 drops_inside = 78462 drops_total = 100000 n = 100000 pi = 3.1384799999999999 Err = 3.1126535897931795E-003 drops_inside = 3674342 drops_total = 4678343 n = 4678343 pi = 3.1415755535667222 Err = 1.7100023070870662E-005`

Other thing is interesting too: With single call

## CODE

## CODE

## RE: Calculating pi

## CODE

and got the following results:

## CODE

## RE: Calculating pi

In my opinion the fact that with n = 4678343 we get a better result than with n = 2147483647 could be possibly due to the fact that also in this plot http://upload.wikimedia.org/wikipedia/commons/7/7d... the error drops to zero for x=5, 9, 12, ..., 130 but at x=12 the error immediately grows up again because of the few points used by the program; while at x=130 the error is low also before and after x=130 because the program made lots of steps to reach that precision.

So maybe n = 4678343 could be the case of x=12, while n = 2147483647 is like x=130. But that's just a hypotesis!

The issue with the single or multiple call...I can't explain it at all!

Your last modified code is very interesting, really!

It seems that attempts=8e13*accuracy, so it looks like they're directly proportional. Nice!

## RE: Calculating pi

## RE: Calculating pi

Monte-carlo methods are based on probabilities ... You may be above or below the limit value after any number of trials with an equal probability. Don't be surprised if you get a better result with a lower number of trials... This is "normal" because the attempts will roughly oscillate around the limit and therefore, from time to time, coming very close to the limit ... by chance !

François Jacq

## RE: Calculating pi

Thanks for the explanation, I was faced with such a problem first time since I've never used Monte Carlo before.