Calculating pi

[pioTx]

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.

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

 

[pioTx]

sorry, I meant an error threshold of 10^(-10)

 

[mikrom]

I found this article

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

and tried this

[COLOR=#a020f0]program[/color] pi_montecarlo
  [COLOR=#2e8b57][b]implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]
  [COLOR=#2e8b57][b]double precision[/b][/color], [COLOR=#2e8b57][b]parameter[/b][/color]::pi[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]3.141592653589793238462643383279[/color]
  [COLOR=#2e8b57][b]double precision[/b][/color] :: pi_approx, err
  [COLOR=#2e8b57][b]integer[/b][/color] :: n
  
  n [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]100[/color]
  [COLOR=#008080]call[/color] compute_pi

  n [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]1000[/color]
  [COLOR=#008080]call[/color] compute_pi

  n [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]100000[/color]
  [COLOR=#008080]call[/color] compute_pi

  n [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]10000000[/color]
  [COLOR=#008080]call[/color] compute_pi
  
  n [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]1000000000[/color]
  [COLOR=#008080]call[/color] compute_pi

  n [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]2147483647[/color]
  [COLOR=#008080]call[/color] compute_pi
  
[COLOR=#a020f0]contains[/color] 

[COLOR=#a020f0]subroutine[/color] compute_pi
  pi_approx [COLOR=#804040][b]=[/b][/color] pi_monte_carlo(n)
  err [COLOR=#804040][b]=[/b][/color] [COLOR=#008080]abs[/color](pi [COLOR=#804040][b]-[/b][/color] pi_approx)
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'n = '[/color], n, [COLOR=#ff00ff]'pi = '[/color], pi_approx, [COLOR=#ff00ff]'Err = '[/color], err
[COLOR=#a020f0]end subroutine[/color] compute_pi

[COLOR=#2e8b57][b]double precision[/b][/color] [COLOR=#a020f0]function[/color] pi_monte_carlo(total_drops)
  [COLOR=#2e8b57][b]implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]
  [COLOR=#2e8b57][b]integer[/b][/color], [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: total_drops
  [COLOR=#2e8b57][b]integer[/b][/color] :: drops_inside, nr_drops  
  [COLOR=#2e8b57][b]double precision[/b][/color] :: x, y

  drops_inside [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color]
  nr_drops [COLOR=#804040][b]=[/b][/color] total_drops
  [COLOR=#804040][b]do[/b][/color] [COLOR=#804040][b]while[/b][/color] (nr_drops [COLOR=#804040][b]>[/b][/color] [COLOR=#ff00ff]0[/color])
    x[COLOR=#804040][b]=[/b][/color]rand()
    y[COLOR=#804040][b]=[/b][/color]rand()
    [COLOR=#0000ff]!write(*,*) 'x =', x, 'y =', y[/color]
    [COLOR=#804040][b]if[/b][/color] ((x[COLOR=#804040][b]*[/b][/color]x [COLOR=#804040][b]+[/b][/color] y[COLOR=#804040][b]*[/b][/color]y) [COLOR=#804040][b]<=[/b][/color][COLOR=#ff00ff]1[/color].) [COLOR=#804040][b]then[/b][/color]
       drops_inside [COLOR=#804040][b]=[/b][/color] drops_inside [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]  
    [COLOR=#804040][b]end if[/b][/color]
    nr_drops [COLOR=#804040][b]=[/b][/color] nr_drops [COLOR=#804040][b]-[/b][/color] [COLOR=#ff00ff]1[/color]
  [COLOR=#804040][b]end do[/b][/color] 

  [COLOR=#0000ff]!write(*,*) 'drops_inside = ', drops_inside[/color]
  [COLOR=#0000ff]!write(*,*) 'total_drops = ', total_drops [/color]

  [COLOR=#0000ff]! return result[/color]
  pi_monte_carlo [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]4[/color]. [COLOR=#804040][b]*[/b][/color] drops_inside [COLOR=#804040][b]/[/b][/color] total_drops
[COLOR=#a020f0]end function[/color] pi_monte_carlo
[COLOR=#a020f0]end program[/color]

Output:

n =        100  pi =  2.8399999999999999   Err =  0.30159274101257338
n =       1000  pi =  3.2120000000000002   Err =  7.0407258987426946E-002
n =     100000  pi =  3.1383999999999999   Err =  3.1927410125733857E-003
n =   10000000  pi =  3.1418368000000001   Err =  2.4405898742685395E-004
n = 1000000000  pi =  3.1415918600000001   Err =  8.8101257311734571E-007
n = 2147483647  pi =  3.1415878809716498   Err =  4.8600409234822450E-006

 

[FJacq]

@pioTx : all your real numbers should be double precision values. Even your declaration

double precision, parameter:i=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:i=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

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

 

[pioTx]

Whoa thank you very much!
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.

 

[mikrom]

Another problem seems to be that n = 2147483647 is 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.

 

[pioTx]

@mikrom: and how could you do that? In fact, in another situation as well I needed a number bigger than 2^31-1, but I didn't manage to solve it. Is there any way?

 

[mikrom]

for example like this

http://gcc.gnu.org/onlinedocs/gfortran/SELECTED_005fINT_005fKIND.html

 

[pioTx]

Ohh thank you! I tried and it worked!

 

[mikrom]

Instead of

x=rand()
y=rand()

I tried to use

call random_number(x)
call random_number(y)

and got these results

$ pi_montecarlo
n =        100  pi =  3.1200000000000001  Err =  2.1592741012573136E-002
n =       1000  pi =  3.0920000000000001  Err =  4.9592741012573160E-002
n =     100000  pi =  3.1394799999999998  Err =  2.1127410125734158E-003
n =   10000000  pi =  3.1424259999999999  Err =  8.3325898742669935E-004
n = 1000000000  pi =  3.1415850000000001  Err =  7.7410125731702806E-006
n = 2147483647  pi =  3.1416497766699876  Err =  5.7035657414328256E-005

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%)

 

[mikrom]

And @pioTx thanks for interesting question:
Although I have heard about Monte Carlo method, but until now I have never used it.

 

[FJacq]

As usual in Monte Carlo method, the precision is proportional to the squared root of n.

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

François Jacq

 

[pioTx]

@mikrom: you're right! call random_number works better than rand(), even though the convergence is very slow.
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/Errore_montecarlo_calcolo_pi_greco.PNG

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

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

and got this

steps=     4678342 pi=   3.1415926536382646      error=   4.8471449076714634E-011

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.

 

[mikrom]

But I wonder why your result with n = 4678343 is more better then my with n = 2147483647:

n =    4678343  pi =  3.1415926536382646   Err =  4.8471449076714634E-011
n = 2147483647  pi =  3.1416497766699876   Err =  5.7035657414328256E-005

With n = 4678343 my code only computes this poor approximation

n =    4678342  pi =  3.1420182620253074   Err =  4.2552101273418685E-004

Is there any error in my code?

 

[mikrom]

If I use rand() instead of random_number() and call the function only once i.e.:

n = 4678343
call compute_pi

I get the desired result:

drops_inside =  3674362
drops_total  =  4678343
n = 4678343  pi = 3.1415926536382646  Err = 4.8471449076714634E-011

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

n = 1000
call compute_pi

n = 100000
call compute_pi

n = 4678343
call compute_pi

I get worse result for n = 4678343

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
[highlight]drops_inside = 3674342
drops_total  = 4678343
n = 4678343  pi = 3.1415755535667222  Err =  1.7100023070870662E-005[/highlight]

Other thing is interesting too: With single call

n = 2147483647
call compute_pi

it delivers worse result then with smaller n = 4678343 :

drops_inside =  1686627151
drops_total  =  2147483647
n = 2147483647  pi = 3.1415878828342949  Err = 4.7707554982068245E-006

 

[mikrom]

I modified the code of pioTx to find out how many attempts meet the requested accuracy:

[COLOR=#a020f0]program[/color] pi_montecarlo
[COLOR=#2e8b57][b]implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]
[COLOR=#2e8b57][b]double precision[/b][/color] :: x, y, error, pi_calculated, accuracy
[COLOR=#2e8b57][b]integer[/b][/color]:: k, n, n_max, circle
[COLOR=#2e8b57][b]double precision[/b][/color], [COLOR=#2e8b57][b]parameter[/b][/color] :: pi [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]3.141592653589793d0[/color]

k [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color]
circle [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color]
n_max [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]2147483647[/color]
accuracy [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]1e-10[/color]
[COLOR=#804040][b]do[/b][/color] n [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]1[/color], n_max
  x[COLOR=#804040][b]=[/b][/color]rand()
  y[COLOR=#804040][b]=[/b][/color]rand()

  [COLOR=#804040][b]if[/b][/color] (x[COLOR=#804040][b]*[/b][/color]x [COLOR=#804040][b]+[/b][/color] y[COLOR=#804040][b]*[/b][/color]y [COLOR=#804040][b]<=[/b][/color] [COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]then[/b][/color]
    circle [COLOR=#804040][b]=[/b][/color] circle [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]
  [COLOR=#804040][b]end if[/b][/color]

  pi_calculated [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]4.0[/color] [COLOR=#804040][b]*[/b][/color] circle [COLOR=#804040][b]/[/b][/color] n

  error [COLOR=#804040][b]=[/b][/color] [COLOR=#008080]abs[/color](pi_calculated [COLOR=#804040][b]-[/b][/color] pi)

  [COLOR=#804040][b]if[/b][/color] (error [COLOR=#804040][b]<=[/b][/color] accuracy) [COLOR=#804040][b]then[/b][/color]
    [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'n ='[/color], n, [COLOR=#ff00ff]'pi ='[/color], pi_calculated, [COLOR=#ff00ff]'Err= '[/color], error
    k [COLOR=#804040][b]=[/b][/color] k [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]
  [COLOR=#804040][b]end if[/b][/color]
[COLOR=#804040][b]end do[/b][/color]
[COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Attempts found ='[/color], k
[COLOR=#a020f0]end program[/color]

and got the following results:

accuracy	Attempts found
1e-10		7989
1e-11		 779
1e-12		  81
1e-13		   8

 

[pioTx]

@mikrom: yeah you're right, I noticed that!
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/Errore_montecarlo_calcolo_pi_greco.PNG

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!

 

[mikrom]

In my naive idea I thought that with the increasing N (number of trials), it increases also the probability that the drop falls into the circle. Either my assumption is wrong, or there is a technical problem - for example that the generated random numbers are not random enough. I have no idea .. but the problem is very interesting

 

[FJacq]

@mikrom your assumption is wrong AND the generated random numbers are probably not random enough.

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

 

[mikrom]

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