×
INTELLIGENT WORK FORUMS
FOR COMPUTER PROFESSIONALS

Log In

Come Join Us!

Are you a
Computer / IT professional?
Join Tek-Tips Forums!
  • Talk With Other Members
  • Be Notified Of Responses
    To Your Posts
  • Keyword Search
  • One-Click Access To Your
    Favorite Forums
  • Automated Signatures
    On Your Posts
  • Best Of All, It's Free!
  • Students Click Here

*Tek-Tips's functionality depends on members receiving e-mail. By joining you are opting in to receive e-mail.

Posting Guidelines

Promoting, selling, recruiting, coursework and thesis posting is forbidden.

Students Click Here

Jobs

Calculating pi

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.

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

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

RE: Calculating pi

I found this article
http://de.wikipedia.org/wiki/Kreiszahl#Statistisch...
and tried this

CODE

program pi_montecarlo
  implicit none
  double precision, parameter::pi=3.141592653589793238462643383279
  double precision :: pi_approx, err
  integer :: n
  
  n = 100
  call compute_pi

  n = 1000
  call compute_pi

  n = 100000
  call compute_pi

  n = 10000000
  call compute_pi
  
  n = 1000000000
  call compute_pi

  n = 2147483647
  call compute_pi
  
contains 

subroutine compute_pi
  pi_approx = pi_monte_carlo(n)
  err = abs(pi - pi_approx)
  write(*,*) 'n = ', n, 'pi = ', pi_approx, 'Err = ', err
end subroutine compute_pi

double precision function pi_monte_carlo(total_drops)
  implicit none
  integer, intent(in) :: total_drops
  integer :: drops_inside, nr_drops  
  double precision :: x, y

  drops_inside = 0
  nr_drops = total_drops
  do while (nr_drops > 0)
    x=rand()
    y=rand()
    !write(*,*) 'x =', x, 'y =', y
    if ((x*x + y*y) <=1.) then
       drops_inside = drops_inside + 1  
    end if
    nr_drops = nr_drops - 1
  end do 

  !write(*,*) 'drops_inside = ', drops_inside
  !write(*,*) 'total_drops = ', total_drops 

  ! return result
  pi_monte_carlo = 4. * drops_inside / total_drops
end function pi_monte_carlo
end program 

Output:

CODE

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 

RE: Calculating pi

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

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

(OP)
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.

RE: Calculating pi

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.

RE: Calculating pi

(OP)
@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?

RE: Calculating pi

(OP)
Ohh thank you! I tried and it worked!

RE: Calculating pi

Instead of

CODE

x=rand()
y=rand() 
I tried to use

CODE

call random_number(x)
call random_number(y) 
and got these results

CODE

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

RE: Calculating pi

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

RE: Calculating pi

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

RE: Calculating pi

(OP)
@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...
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

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

CODE -->

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.

RE: Calculating pi

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

CODE

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

CODE

n =    4678342  pi =  3.1420182620253074   Err =  4.2552101273418685E-004 
Is there any error in my code?

RE: Calculating pi

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

CODE

n = 4678343
call compute_pi 
I get the desired result:

CODE

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:

CODE

n = 1000
call compute_pi

n = 100000
call compute_pi

n = 4678343
call compute_pi 
I get worse result for n = 4678343

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

n = 2147483647
call compute_pi 
it delivers worse result then with smaller n = 4678343 :

CODE

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

RE: Calculating pi

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

CODE

program pi_montecarlo
implicit none
double precision :: x, y, error, pi_calculated, accuracy
integer:: k, n, n_max, circle
double precision, parameter :: pi = 3.141592653589793d0

k = 0
circle = 0
n_max = 2147483647
accuracy = 1e-10
do n = 1, n_max
  x=rand()
  y=rand()

  if (x*x + y*y <= 1) then
    circle = circle + 1
  end if

  pi_calculated = 4.0 * circle / n

  error = abs(pi_calculated - pi)

  if (error <= accuracy) then
    write(*,*) 'n =', n, 'pi =', pi_calculated, 'Err= ', error
    k = k + 1
  end if
end do
write(*,*) 'Attempts found =', k
end program 

and got the following results:

CODE

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

RE: Calculating pi

(OP)
@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... 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

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 smile

RE: Calculating pi

@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

RE: Calculating pi

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

Red Flag This Post

Please let us know here why this post is inappropriate. Reasons such as off-topic, duplicates, flames, illegal, vulgar, or students posting their homework.

Red Flag Submitted

Thank you for helping keep Tek-Tips Forums free from inappropriate posts.
The Tek-Tips staff will check this out and take appropriate action.

Reply To This Thread

Posting in the Tek-Tips forums is a member-only feature.

Click Here to join Tek-Tips and talk with other members!

Close Box

Join Tek-Tips® Today!

Join your peers on the Internet's largest technical computer professional community.
It's easy to join and it's free.

Here's Why Members Love Tek-Tips Forums:

Register now while it's still free!

Already a member? Close this window and log in.

Join Us             Close